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Abstract. Core collapse supernovae are a promising source of detectable gravitational waves. Most of the existing 
(multidimensional) numerical simulations of core collapse i n general relativity h ave be en done using approxima- 
tions of the Einstein field equations. As recently shown bv iDimmelmeier et alJ i2002al lbh. one of the most inter- 
esting such approximation is the so-called conformal flatness condition (CFC) of Isenberg, Wilson and Mathews. 
Building on this previous work we present here new results from numerical simulations of relativistic rotational 
core collapse in axisymmetry, aiming at improving the dynamics and the gravitational waveforms. The computer 
code used for these simulations evolves the coupled system of metric and fluid equations using the 3+1 formalism, 
specialized to a new framework for the gravitational field equations which we call CFC+. In this approach we add 
new degrees of freedom to the original CFC equations, which extend them by terms of second post-Newtonian or- 
der. The resulting metric equations are still of elliptic type but the number of equations is significantly augmented 
in comparison to the original CFC appr oach. The hydrodynamics e volution and the CFC spacetime metric are 
calculated using the code developed by IDimmelmeier et alJ J2002a) , which has been conveniently extended to 
account for the additional CFC+ equations. The corrections for CFC+ are computed solving a system of elliptic 
linear equations. The new formalism is assessed with time evolutions of both rotating neutron stars in equilibrium 
and gravitational core collapse of rotating polytropes. Gravitational wave signals for a comprehensive sample of 
collapse models are extracted using either the quadrupole formula or directly from the metric. We discuss our 
results on the dynamics and the gravitational wave emission through a detailed comparison between CFC and 
CFC+ simulations. The main conclusion is that, for the neutron star spacetimes analyzed in the present work, no 
significant differences are found among CFC, CFC+, and full general relativity, which highlights the suitability 
of the former. 
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1. Introduction 

One of the most interesting topics in relativistic astro- 
physics is the formation and evolution of compact objects, 
such as those produced in the rotational core collapse of 
massive stars, or in the merging of a compact binary con- 
sisting of two white dwarfs, two neutron stars, or a white 
dwarf and a neutron star. The catastrophic events in which 
the new-born compact objects are formed are a promising 
source of gravitational waves in the kHz frequency range, 
hence optimal to be detected at galactic distances with the 
new ground-based interferometer detectors such as LIGO 
or VIRGO. 

Send offprint requests to: Pablo Cerda-Duran, 
e-mail: pablo.cerda@uv.es 



Gravitational radiation is produced by accelerating 
matter whose dynamics is not spherically symmetric. 
Among the list of astrophysical sources of gravitational 
radiation, axisymmetric sources such as rotational core 
collapse are particularly amenable to numerical investi- 
gations, as present-day computational resources allow for 
more accurate calculations than in full three dimensions. 
The early phase of the merging of two compact objects, 
when the two bodies start to plunge from the innermost 
stable circular orbit, is a three-dimensional situation, al- 
though after the merging itself the scenario is again suit- 
able for investigations in axisymmetry. On the other hand, 
Newtonian gravity does not describe correctly the dynam- 
ics of these processes, and one has to use general relativity 
for a proper description, as well as to compute waveforms, 
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or at least some approximation better than the Newtonian 
one. Post-Newtonian approximations of the metric in the 
near zone match quite well with the features described 
before, and probably most of these scenarios can still be 
handled without fully general relativistic calculations: Let 
d be the typical dimension of the source, v the typical ve- 
locity of the system, and t the typical timescale of varia- 
tion of the gravitational field; hence, the typical frequency 
of the gravitational waves in the wave zone is 1 /t, and the 
wavelength is A ~ ct, where c is the speed of light. From 
the virial theorem it follows that t is also the typical dy- 
namical timescale, hence d ~ vt and v/c ~ d/ct ~ d/X. 
Therefore, a post-Newtonian expansion around v/c — ► 
is equivalent to an expansion around d/X — > 0. For the 
typical scales involved in core collapse and neutron stars 
we get d 2 /A 2 ~ 0.05. That estimate yields first and second 
order post-Newtonian corrections of about 5% and 0.3%, 
respectively. We note, however, that although these cor- 
rections may seem to be small, the non-linearity of the 
equations makes it difficult to predict beforehand the ef- 
fects in a numerical simulation. 

During the last two decades calculations of gravi- 
tational waveforms from core collapse ha ve been done 
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200fl IShibata fc Sekiguchil l2004at lOtt et alJl2004h . First 
attempts were made in Newtonian g ravity, using 
Eulerian c o des with artificial visco sity 
iFirml Il989t iMonchmever et all Il99ll) to calculate the 
hydrodynamics. Pse u do-sp ectral methods were used by 
iBonazzola fe Marckl ljl993|) . who followed the collapse 
without including the bounce phase due to the appearance 
of numerical instabilities associated with the presence of 
a shock wave. The need for correctly treating the shock 
wave, which forms after core bounce, gradually led to the 
use of high-resolu t ion sh ock-capturing (HRSC) schemes. 
IZwerger fe Mullerl l|l997l) first used HRSC methods to sim- 
ulate a sequence of collapsing rotating polytropes in ax- 
isymmetry, providing a comprehensi ve waveform catalog . 
Extensions to 3D were carried out bv lRampp et al.l (fl9 98"> . 
In bo th works a simple equation of state l|Janka et alJ 
Il993|) was used. Recently there have been attempts 
to inc lude more realist ic physics in Newtonian simula- 
tions: lOtt et al.l ll2004h used t he equation of sta te of 
lLattimer fc Swestvl l|l99ll . while IMiiller et all l|2004h also 
incorporated 2D neutrino transport and computed the 
gravitational wave emission produced by neutrino-driven 
convection. 

A step forward to improve the treatment of th e grav- 
itationa l field was taken by iDimmelmeier et"aTI l|200ll 
l2002alrJ) . who used the conformal flatness condit ion (CFC 
hereaf ter) of Isenberg, W i lson, a nd Mathews (see llsenberd 
l|l978|) and lWilson et alJ |l990) ). Note that in the spher- 
ical case CFC is not an approximation to general relativ- 
ity but rather is exact. In a more general context it has 



been shown that conformally flat gravity agrees with full 
general rel ativity at the level of the first post-Newtonian 
expansion l)Klev fe Schaferlll999j) . Therefore the CFC ap- 
proximation reproduces quite well the main features of 
the models studied (quasi-spherical sources tractable by 
the post-Newtonian formalism). Numerically, the major 
advantage of the CFC approximation is that the 3+1 met- 
ric equations are reduced to elliptic equations, which, in 
principle, makes their implementation more stable than in 
the standard formulations of the field equations commonly 
used in nu merical relativity. Co n cerning the core collapse 
waveforms, IDimmelmeier et"al1 ((2002b) found quantita- 
tive and also qualitative differences between Newtonian 
and CFC calculations. This demonstrates that a relativis- 
tic calculation of the spacetime dynamics is necessary 
for properly studying core collapse and the gravitational 
waves emitted. Recentl y this work in axisymmet ry has 



been extended to 3D bv IDimmelmeier et al. (|200 



com- 
bining HRSC methods for the hydrodynamic equations 
and spectral methods for the metric equations. 

First multidimensional core collap se simulation s in ful l 
general relativity were reported by ISiebel et all l|2003h 
who used the characteristic formulation of the Einstein 
equations. In these simulations the gravitational wave- 
form extraction using the Bondi news function at future 
null infinity was hampered by gauge effects and numer- 
ical inaccuracies. Recently, IShibata'fc Sekiguchil ll2004ah 
have performed simulations in axisymmetry using the 
so-called Cartoon method to solve the Einstein equa- 
tions, fi nding remarkable a g reemen t with the CFC re- 
sults of IDimmelmeier et"al"l ( 2002bl ) . Extensions of this 
work to full 3D hav e been reported inlShibata fc Sekiguchil 
l)2004b|) . Except for lSiebel etafl l)2003|) . in all existing sim- 
ulations thus far the gravitational wave signals have been 
computed using the Newtonian quadrupolc formula. 

In this paper we propose a new approximation scheme 
for the metric equations that is an extension of the CFC 
approach. Our approximation consists of adding second 
post-Newtonian terms to the conformal flet metric. With 
the new approach, which we call CFC+, the resulting met- 
ric equations are still elliptic, and reduce to full general 
relativity in the spherical case. Therefore, we conserve the 
main advantages of the CFC approximation at small extra 
computational cost (the additional elliptic equations are 
linear), while improving the computation of the dynamics 
and the metric. The aim of this paper is to compare the 
quality of the CFC+ approximation in different scenarios 
with that of the CFC approach. In particular we focus on 
the study of two scenarios of current interest in relativistic 
astrophysics, the gravitational collapse of rotating stellar 
cores to neutron stars (the relevant mechanism responsi- 
ble for core collapse supernovae of type II/Ib/Ic), and the 
stability and oscillations of rotating neutron stars. 

This paper is organized as follows: In Sect. |21 we 
present the mathematical framework and introduce the 
general relativistic equations for the hydrodynamics and 
the spacetime metric that we implement in the numerical 
code. Particular emphasis is given to the derivation of the 
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CFC+ equations (for technical details see Appendices 1X1 
IbI andEJ). Sect. |21 is devoted to describing the numeri- 
cal methods used to solve the hydrodynamics and metric 
equations. We further explain in this section the proce- 
dure we follow to extract gravitational waves (see also 
Appendix lD|) . In Sect.EJwe specify the models used for the 
simulations of rotating relativistic stars. The results from 
the simulations of evolutions of rotating neutron stars and 
rotational stellar core collapse to a neutron star are dis- 
cussed in Sects. El and El respectively. The conclusions of 
this work are summarized in Sect.0 

We use a spacelike signature (— , +, +, +) and units in 
which c = G = 1 (except in some passages, particularly in 
the Appendices, where c and G are retained for a better 
understanding of the post-Newtonian expansion). Greek 
indices run from to 3, Latin indices from 1 to 3, and 
we adopt the standard Einstein convention for the sum- 
mation over repeated indices. The indices of 3- vectors and 
3-tcnsors are raised and lowered by means of the 3-mctric. 

2. Mathematical framework 

We adopt the 3 + 1 formalism (jArnowitt et alJ Il962|) 
to foliate spacetime with a metric g^ into a set of 
non-intersecting spacelike hypersurfaces. The line element 
reads 



ds 2 



-a 2 dt 2 + jijidx* + ^dt)(dx j + f3 j dt), 



(1) 



with a being the lapse function, which describes the rate of 
advance of time along a timelike unit vector normal to 
a hypersurface, /3 1 being the spacelike shift 3- vector, which 
describes the motion of coordinates within a surface, and 
7y being the spatial 3-metric. 

2.1. Hydrodynamic equations 

The hydrodynamic evolution of a relativistic perfect fluid 
with rest-mass current J M = pu 11 and energy-momentum 
tensor T^ v — phu^u u + Pg^ v in a (dynamic) spacetime is 
determined by a system of local conservation equations, 
which read 



V„ J" = 0, VuT^ = 0, 



(2) 



where the covariant derivative is denoted by V^. The fluid 
is specified by the rest mass density p, the specific enthalpy 
h = 1+e+P/p, the pressure P, the specific internal energy 
e, and the 4- velocity m m . The pressure is determined by an 
equation of state (EOS) P = P(p, e). 

The 3-velocity of the fluid, as measured by an Eulerian 
observer at rest in a spacelike hypersurface S t , is 



+ 



(3) 



Following the work of iBanvuls et alJ l|l997|) we introduce 
a set of variables in terms of the primitive (physical) hy- 
drodynamic variables (p, v l , e) , whose associated densities 



are conserved: 

D = pW, 
S l = P hW 2 v\ 
t = phW 2 - P - D. 



(4) 
(5) 
(6) 



In the above expressions W is the Lorentz factor defined as 
W — au* , which satisfies the relation W = 1/ Vl — ViV 1 . 

The previous choice allows us to write the conserva- 
tion laws J2Jl as a fir st-order, flux-conserv ative hyperbolic 
system of equations l|Banvuls et al.lll99^) : 



1 



-9 



dt 



dx l 



= Q, 



(7) 



with the state vector, flux vector, and source vector given 

by 



U=[D,S j ,r], 

F' 1 = [Dv\ SjV 1 + 8iP, TV 1 + Pv l ] 

, d In a 



Q 



'2 dxi ' 



rjifJ,0 



rpfAV pi) 



(8) 
(9) 

(10) 



Here v % = v % — f3 l /a, and \/—g = a<y/7, with g = det(<7 A1!/ ) 
and 7 = det(7y) being the determinant of the 4-metric 
and 3-metric, respectively; r* v are the Christoffel sym- 
bols. 

For convenience, we also define the following modified 
conserved quantities: 



D* 



VrD, s* = yfiSi 



(ii) 



The determinant 7 is actually the ratio of 7 evaluated 
on some coordinate grid {x^} and the determinant 7 of 
the flat metric 7^ on the same grid: 7 = 7/7. 

2.2. 3+1 metric equations 

In the 3 + 1 formalism, the Einstein equations split into (i) 
evolution equations for the 3-metric 7^ and the extrinsic 
curvature Ky, 

dan = -2aK tJ + 2V(i/3j) , (12) 
d t K l3 = -ViV.Q + a(Rij + KKij - 2K ik K*) 
+ (] k V k K l] +2K k{t V 3) l3 k 

(13) 



Swafa-ljfiSZ-E!) 



and (ii) constraint equations, 

R + K 2 - /\,,/\'' - 16ttE = 0, 
Vi(K ij ~Y J K)~8TrS j = 0, 



(14) 
(15) 



which must be fulfilled at every spacelike hypersurface. 

In these equations Vj is the covariant derivative with 
respect to the 3-metric 7^, K+j is the extrinsic curvature, 
Rij is the Ricci tensor on E t , R the scalar curvature and 
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K = K\ the trace of the extrinsic curvature. The brack- 
ets around indices indicate symmetrization. The projec- 
tion of the energy-momentum tensor onto the spatial hy- 
persurface is defined as Sij = phW 2 ViVj + lijP, while 
E = phW 2 — P is the total energy per volume unit as 
measured by an observer at rest in E t . 

The ADM gauge, which we will be using in the deriva- 
tion of the metric equations in Sects. I2~3l and l2.4l is defined 
as that gauge for which the 3-metric can be decomposed 
into a conformally flat term plus a transverse traceless 
term, 



lij = 
with 



J.4 - 
<P lij 



. TT 



0, 



? fc v fc fcF = o, 



(17) 



= 

where 7™ is the flat metric (with inverse j tJ ; Ty = ^ = 
in Cartesian coordinates), 4> is the conformal factor, 
and hj^ is transverse and traceless. The conjugate mo- 
mentum tt^ of 7y- is traceless as well: 

7T lJ %- = 0. 



(18) 



2.3. CFC metric equations 



In spherical symmetry hj^ = 0, i.e. the 3-metric is con- 
formally flat (CF). A first reasonable approximation in 
quasi-spherical scenarios is imposing a vanishing hj^ in 
Eq. JEi): 



lij = lij- 



(19) 



This is the conformal flatness condi tion or I se nberg- 
Wilson-Mathews approximation ijlsenberd [l978; 
IWilson et all Il996j) . In the ADM gauge a confor- 
mally flat 3-metric ensures the maximal slicing condition, 
K = 0. Hence, it can be shown that the 3 + 1 metric 
equations (|12H15J1 reduce to a set of five coupled elliptic 
(Poisson-like) equations for the metric components, 



K K lJ 

Acj) = -2^ I phW 2 -P+ 

' ' 1D7T 



(20) 



A{a<£) = 2ira<j) 5 [ ph(3W 2 - 2) + 5P + , (21) 

\ 1D7T J 

A/3 1 = Wira^S 1 + 2cf) W K i: >ij ] (^] ~ iv l V fe /3 fe , (22) 



where V and A = 7 y V,Vj are the flat space Nabla and 
Laplace operators, respectively. 

Applying the CFC to the traceless part of Eq. (jT2")l 
yields the following relation between the conformal factor 
and the shift vector components: 



d t 4> = f V fc/ 3 fc . 
6 



(23) 



With this, Eq. (|12|l transforms into an expression for the 
extrinsic curvature components: 



Kij = (yiPj + - V fe /3 fe 



(24) 



As a consequence of setting the off-diagonal elements 
of 7^ to zero, the degrees of freedom representing gravita- 
tional waves are removed from the spacetime. Therefore, 
we calculate gravitational waveforms in a post-processing 
step using the quadrupole formula. Note that at that ap- 
proximation level the ADM gauge is equivalent to the 
maximal slicing quasi-isotropic (MSQI) gauge, i.e. max- 
imal slicing (K — 0) plus quasi-isotropic coordinates 
(g r g = and geg = g rr in spherical coordinates in an 
orthonormal basis). 



(16) 2.4. CFC+ metric equations 



We improve the metric approximation by adding to the 
CFC metric the second post-Newtonian deviation from 
isotropy. We call this approximation CFC+. Compared to 
the CFC metric (O the CFC+ metric 



lij 



lij 



,2PN 



(25) 



now includes a new traceless and transverse term /i^ PN , 
which is identical to hj^ in the decomposition 1)16(1 up 
to the second post-Newtonian expansion. As described in 
Appendix^] this higher order correc tion h 2 f^ is t he solu- 
tion of the tensor Poisson equation llSchafe'rlll990Tl 



Aft H =lii tkh 



where the source F^i is given by 

St St 



Fm = ~4\7 k UX7,U- 16ir 



D* 



(26) 



(27) 



with U being the "Newtonian" potential, which is the so- 
lution of the Poisson equation 



AU = -AttD* 



(28) 



The operator 7^™ is the transverse, trace-free (TT) pro- 



jector as defined in Eq. 1IA.6I) . 

Of the original CFC metric equations H2UI - I22|) for cf>, a, 
and (3 l , only the equation for the lapse function, Eq. (|21|) , 
has to be modified. It now depends explicitly on the second 
post-Newtonian correction /if™: 



A(a<f>) = A (a 



1 1 K 



(29) 



The metric equations for the conformal factor (f> and the 
shift vector /3* remain unaltered. However, the CFC+ cor- 
rections couple implicitly as the source terms on their 
right-hand sides depend on a. 

The computation of /if™, i.e. the inversion of the ten- 
sor Poisson equation (|26() by means of the Poisson integral 
operator A , is simplified by the introduction of interme- 
diate potentials S, Si,Ti,lZi, and Sij, which are solutions 
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of the following scalar/vector/tensor-Poisson equations: 

(30) 



AS = -47r^-^-^a: J , 



A5, = 
A7^ = 



-4-7T- 



D* 
S*S* 

D* 



ATZi = V l {V 1 UV k Ux 3 x k ), 

s*s* 



ASij = 



-47T- 



ViUVjU. 



(31) 

(32) 
(33) 
(34) 



These equations are designed such that their source terms 
approach zero like r -3 when r = \x\ tends towards in- 
finity, which ensures the existence of the corresponding 
Poisson integrals. With the help of these potentials h 2 J n 
can simply be expressed as 



2PN 



= A" 



("TTM t? \ 



3x k V {i S j)k + -7 Jm a; m V 4 (j kl S kl ) 



1 



-.i J Y,Y.Si, 



■la 



\l kl S k i + x k f m V m S kl - l U ^kSi 



(35) 



as shown in Appendix iBl 

Due to the specific type of elliptic solvers employed in 
our computer code (see Sect. 13.21 below), it is not possi- 
ble to use the inverse image method to evaluate the in- 
termediate potentials up to spatial infinity. We instead 
solve Eqs. (I3()H34(I assuming specific boundary conditions. 
These are determined from the multipole expansion M. 
of the intermediate potentials S, Si, %, IZi, and Sij. As 
detailed in Appendix El these multipole expansions are 
given by 

M (S) = - f d 3 x ^ ( ^x k x) , (36) 



M(Si) 



M(%) 



r 
1 

r 

A£2 

~2~r~ 
1 

r 

M 2 
~2r 



d 3 x V7D 



D 

SfS 



D* 



d 3 x 



k x k +%x j (U + x k X7 k U) 
(37) 



7 £T 



l kl S* k Sf 

D *2 



+ U)%3t l 



M(Ki) 



M 2 ^n 3 



(38) 
(39) 

+ \%U + % k x k djU\, 
(40) 

modulo 0(l/r 2 ) corrections. M = J d 3 Xy/^D* is the 
baryonic rest mass. By evaluating the above expansions 
at a finite radius outside the star we obtain the necessary 
boundary conditions for S, Si, %, IZi, and Sij. 



MiSfj) = i Jd 3 x^D* 



sjs; 

D* 2 



3. Numerical implementation 

The code used for the simulations pre sented in this paper 
is an exte nsion of the one described in lDimmelmeier et alJ 
( 2002a b), incorporating the new CFC+ metric equations. 
Hence, the interested reader is addressed to those refer- 
ences for details additional to those discussed below. 

We use spherical polar coordinates {t, r, 9, (f} and as- 
sume axial symmetry with respect to the rotation axis and 
symmetry with respect to the equatorial plane at 9 = w/2. 
The computational grid is composed of n r radial zones and 
ng equidistant angular zones, whose specific values for the 
simulations reported here are given below. For convenience 
the radial cell-spacing is chosen equidistant for evolutions 
of equilibrium neutron stars and logarithmically spaced 
when simulating core colla pse. As in the original code 
l|Dimmelmeier et alJ I20021H) the part of the grid outside 
the star is filled with an artificial atmosphere (as custom- 
ary in finite difference cod e s similar to ours; seelFont et alJ 
l|2002l) : lDuez et alJ l|2002l) : iBaiotti et all \2Q(\£\ ). This at- 
mosphere obeys a polytropic equation of state and has a 
very low density such tha t its presence does not affect the 
dynamics of the star (see iDimmelmeier et"aTI ()2002al) for 
details). Moreover, an extended grid containing no mat- 
ter is used beyond the atmosphere for the CFC+ metric 
calculations, namely to properly capture the radial fall-off 
behavior of the metric potentials and to extract the grav- 
itational radiation using components of hj^ in the wave 
zone (see Sect. I6.2)l . 



3.1. Hydrodynamics solver 

The hydrodynamics solver performs the numerical inte- 
gration of the system of hydrodynamic conservation equa- 
tions using a high-resolution shock-capturing (HRSC) 
scheme. This method ensures numerical conservation of 
physically conserved quantities and a correct treatment 
of discontinuities such as shocks. In HRSC methods a 
Riemann problem has to be solved at each cell inter- 
face, which requires the reconstruction of the primitive 
variables (p,v l ,e) at this interfaces. We use the PPM 
method of Colella k Woodward! l|l984f) for the reconstruc- 
tion, which yields third order accuracy in space. The so- 
lution of the Riemann problems then provides the numer- 
ical fluxes at cell interfaces. To obtain this solution, the 
characteristic struc ture of the hydrodyn amics equations is 
explicitely needed ( Banvuls et alJll997l) . 

Contrary to CFC in the CFC+ approximation the met- 
ric does not have any zero off-diagonal elements, and thus 
we have to use the most general eigenvalu es and eigen - 
vectors in general relativity as reported in iFontl l|2003|) . 
Once the spectral information is known, the numerical 
fluxes are co mputed by means of Marqu ina's approximate 
flux formula l|Donat k Marauinal ll996). The time update 
of the conserved vector U is done using the method of 
lines in combination with a Rungc Kutta method with 
second order accuracy in time. Once the conserved quanti- 
ties (D, Si,r) are updated in time, the primitive variables 
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need to be recovered. This is done through an iterative 
Newton-Raphson method, as these variables can not be 
obtained in closed form from the conserved variables. 

We note that the sources Q of the hydrodynamic equa- 
tions have been implemented in the code using a general 
form for the metric, although they can be simplified for a 
metric with vanishing nondiagonal terms in the 3-metric, 
as, for example, in CFC. 

3.2. Metric solver 

The main feature of the approximations we are using for 
the metric is that only elliptic equations have to be solved 
to know the metric at each time step. In our approximate 
scheme we solve the equations hierarchically. First a solu- 
tion of the Poisson equation 1(28(1 for U is obtained. Then 
we solve the equations ((30H34J) for the intermediate po- 
tentials S, Si, %, TZi, and Sij, which we need to calculate 
hj^- in Eq. . Finally, we use hj^ to solve the modified 
equations for (j>, (3 l , and a, Eqs. (J2J3 E3 EHl ■ For each step 
we use different methods according to the mathematical 
characteristics of each equation. 

The equation for the Newtonian potential U is a lin- 
ear Poisson equation with compact support sources, and 
hence standard methods for Poisson equations like inte- 
gral solvers can be used to obtain a numerical solution. 
For an equation of the form 



Au(x) = f(x), 



(41) 



the solution for the potential u can be expanded in ax- 
isymmetry as 

«(*) = --!- /dV^A 

Att \x — x'\ 



1 oo 

1=0 

where 



(42) 



(43) 



|x'|<fl 



= r l I dr' d// f&)-j^Pi(tM'). (44) 

J\m'\>R r 

Here Pi are the Legendre polynomials, and fi = cos 9. 
We numerically integrate Eqs. 1(431 144(1 by assuming f(x') 
to be constant inside each computational cell, integrat- 
ing over r' and 9' analytically within each cell, and then 
adding up the partial integrals to obtain the solution at ev- 
ery grid point of the computational d omain. This method- 
has b een d escribed a nd te sted in iMiiller fe Steinmet d 
lll995l) and IZwergerl lll995l) . and successfully used by 
IZwereer fe Miillerl l)l997|) to calculate the Newtonian po- 
tential in axisymmetric core collapse simulations. 

The equations (|30H34|I for the intermediate potentials, 
which are needed to compute the second post-Newtonian 
corrections to CFC, can in general be solved separately 
as a scalar Poisson equation for S, three vector Poisson 



1 



equations for Si,%, and TZi, and a tensor Poisson equation 
for S^. In the axisymmetric case we can take advantage 
of some additional simplifications: The 93-component of 
the vector-Poisson equations decouples from the r- and 
^-components, and the rip- and ^-components of the 
tensor-Poisson equation decouple from all other compo- 
nents, even though e.g. A<S, 7^ (AS)j (similar for the other 
vectors and tensor), which means that the various compo- 
nents in general couple to each other. Therefore, the equa- 
tions can be grouped into 9 sets of linear elliptic equations: 
four sets of one equation, four sets of two equations, and 
one set of four equations, with coupling only within the 
respective set. 

The discretization of each of the equations on the 
{r, (9}-grid leads to 9 sparse linear matrix equations of the 
type 



Au = f, 



(45) 



where u — is the vector of unknowns with i , j la- 
beling the grid points and k ranging from 1 to 1, 2, or 
4 depending on the number of coupled components. The 
vector of sources is respectively denoted as / = /£, and 
^4 is a matrix which does not depend on u, as the origi- 
nal system is linear. The linearity of the equations allows 
us to avoid an iteration scheme and to use instead the 
LU decomposition method to invert A.. The main advan- 
tage of the LU method is that the decomposition itself 
(which is the most time consuming step) only has to be 
done once at the beginning. Then it can be used to calcu- 
late the solution for different source vectors / during the 
metric computations, which is computationally very fast 
and efficient. The LU decomposition is performed using 
standard LAPACK libraries for banded matrices. We use 
the monopole solution of Eqs. ((361 - 140(1 as explicit outer 
boundary condition for the intermediate potentials. This 
procedure is only accurate far from the sources, and the 
matching of the metric to a monopole solution deterio- 
rates if it is performed too close to a strongly gravitating 
nonspherical star. As a consequence when calculating the 
components of the spacetime metric we use a grid which 
extends well beyond both the boundary of the star and 
also the atmosphere, and apply the boundary condition 
at the outer boundary of this extended vacuum grid. An 
example of the influence of the location where the bound- 
ary condition is applied on the accuracy of the metric 
solution is presented in Sect. 15.11 Although the nonlin- 
ear metric equations for </>, a, and (3 l have source terms 
with noncompact support, the location of the outer grid 
boundary has much les s influence on the a c curac y of the 
numerical solution (see iDimmelmeier et al.l 1)20041) ). 

At this stage we have a numerical solution for hj^ 
and are ready to solve the CFC+ metric equations for 
the conformal factor <fi, the shift vector (3 l , and the lapse 
function a, Eqs. 1(2011221129")) . For the comparison between 
CFC and CFC+ presented in Sects. and EJ we also need 
to solve the CFC metric equations, which, we recall, are 
equivalent to the CFC+ equations up to corrections in the 
metric equation for a (|2"TJ) . As both systems of equations 
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arc 5 nonlinear elliptic coupled Poisson-likc equations, we 
can apply the same methods to solve them. In either case, 
we can write them in compact form as 



Au(x) = f(x;u(x)), 



(46) 



where u = u k = ((/>, ct(f>, ft ) , and / = f k is the vector of 
the respective sources. These five scalar equations for each 
component of u couple to each other via the source terms 
which in general depend on the various components of u. 

We use a fix-point iteration scheme in combination 
with the linear Poisson solver (|42|l described above to solve 
these equations. The value of u at each iteration s, de- 
noted by it s , is set constant in the sources / to calculate 
a new value it s+1 , 



x) = A- 1 f(x;u s (x)) 



(47) 



As a result the 5 previously coupled nonlinear equations 
reduce to a decoupled set of 5 linear scalar Poisson equa- 
tions. The solution vector u s+1 at each iteration step is 
obtained by solving the associated 5 Poisson equations of 
the type f4"T|l separately. After the computation of u s+1 
the right-hand side of Eq. I|47[) is updated by replacing 
u s — ► u s+1 , which is used as a new starting value for the 
next iteration. Convergence to the desired numerical solu- 
tion u is achieved when the relative variation \u s+1 /u s — 1| 
of the numerical solution of u between two successive iter- 
ations is smaller than a certain threshold, which we set to 
10 -6 . In the simulations reported here the metric solver 
successfully converges when using the flat metric as initial 
guess in each metric computation. However, to acceler- 
ate convergence, during the hydrodynamic evolution we 
take the metric from the previous metric computation as 
starting value for the subsequent one. Furthermore, we 
cut the Legendre polynomial expansion in Eq. I|42|) at 
I = 10. For the CFC metric equations in axisymmetry 
th e above method ha s been recently discussed in detail 
in iDimmelmeier et all l)2004l) . 

3.3. Gravitational wave extraction 

For an isolated source of gravitational radiation the metric 
can be considered as asymptotically flat. In the wave zone, 
defined as a neighborhood of future radiative infinity X + 
for which r>A/ (2tt), the metric can be decomposed into 
the Minkowski metric r]^ v plus a small linear perturbation 
h/jf?. By an appropriate choice of gauge it is always pos- 
sible to make the latter quantity algebraically transverse 
and traceless, up to quadratic corrections. At the linear 
order we have then 



. rad _ p 



■TTfcZ 



ilkl — 7jm)i 



(48) 



where P?™ = %(p%)j [(nPn k - jP k ) {n q n l - 
(n p n q - Y q ) (n k n l - ■y kl ) /2] stands for the algebraic 
transverse traceless projector. We note that the above 
expression stays invariant with respect to asymptotically 



Minkowskian gauge transformations of the 3-metric. It in- 
volves only two independent components, h + and h x , rep- 
resenting the two degrees of freedom of the gravitational 
waves. 

At linear order in the wave zone, h™ A can be approxi- 
mated by means of the Newtonian quadrupole formula in 
some Bondi-like coordinate system. In this representation, 
the waveform depends only on the second time derivative 
(denoted by a double dot) of the mass quadrupole moment 
Qij of the source: 



^ ad ( 



x ^ _ pTTkl 2Qkl (tret) 



(49) 



with the retarded time i rc t = t — r/c. The argument of 
Qki in the ADM gauge is not simply eq ual to t — r/c . 
It also contains a logarithmic correction llSchafeilll99rl 
-2 (M + C(l/c 2 )) Gln[r/(crj)]/c 3 + 0(G 2 ), where r) de- 
notes an arbitrary constant with the dimension of time 
that cancels a term entering the first tail contribution. For 
large distances near future null infinity X + , this correction 
cannot be neglected, but it may be omitted in investiga- 
tions focusing on time dep endences as it depends on r 
only. Following iFinnl l|1989|) . the time derivatives in the 
standard quadrupole formula can be replaced by spatial 
derivatives making use of the equations of motion. This 
yields the so-called stress formula, 



Qij = 2ld 3 x 



lk<ilj>iv k v l + x k % <i V j> U 



(50) 



where the brackets <> denote the symmetric and trace- 
free sets of indices, e.g. 7fe<iVj>[/ = (jki^jU + 
% j V i U)/2 - % V fc (7/3. 

In axisymmetry, h\^ A has only one degree of freedom, 
h+ . In the quadrupole approximation and in spherical co- 
ordinates, it takes the form 



1 / 1 ^ 

hf M {x,t) = -\ — 



8V 7T 



sin 



A%(t-r/c) 



(51) 



with (x Qgg being the amplitude of the quadrupole 
signal. The above formula allows us to calculate waveforms 
directly from the sources, with no need of knowing the 
dynamical part of the metric. It may be used with profit in 
the CFC approximation, for which the metric is isotropic 
and has therefore no gravitational wave content. 

In the CFC+ approach, the non-conformal part of the 
3-metric in the near zone at the dominant (second) post- 
Newtonian order, /if™, is the solution of the pseudo- 
Poisson equation l|26ll . where the remainder C(l/c 6 ) has 
been removed. This equation has an elliptic character, so 
that the information it carries propagates instantaneously. 
As a result, /i? PN is neither wave-like nor associated to 
any retardation effects, reflecting the fact that the post- 
Newtonian expansion is only valid at distances comparable 
to the wavelength of the system A. By contrast, we know 
that the full non-conformally flat part of the 3-metric 
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carries the dynamical degrees of freedom of the gravita- 
tional field. This indicates a relation between M" n and 
h\j d near I + (they are indeed equal). Now, remarkably, it 
turns out that there is also a link between /t 2PN and the 
radiative metric. As explained in Appendix [Dj the alge- 
braic transverse traceless projection of the former quan- 
tity evaluated at retarded time t ie t behaves as ft.^ uad /8 at 
future null infinity, P™h 2 ™{x,t rct ) ~ ^ uad /8, or in 
components, 



rotating matter distributions, whose pressure obeys the 
polytropic relation 



/>fV,i rct )~i/4 uad 0M), 

hf N ( X ,t Ict )~lhr d (x,t), 



(52) 
(53) 



where /i+ PN can be trivially calculated from the combina- 
tion h+ — (hJJ — ti£p)/2, whereas fo 2PN is equal to h 2 ™ . 

We see that the two polarizations extracted directly 
from hW N in the wave zone of the true waves differ from 
the quadrupole waveforms by (i) a factor 8, and (ii) the ab- 
sence of retardation. The gravitational waveforms can thus 
be evaluated directly from /i+ PN and /i 2PN even though 
there are no propagating gravitational waves in the CFC+ 
spacetime (which just means that there is no radiation 
back-reaction due to energy losses caused by gravitational 
waves). Furthermore, the wave amplitude can be deduced 
solely from /i^_ PN in the axisymmetric case. 

In our simulations we extract the waveforms at a dis- 
tance r ~ 50A/(27r) in the equatorial plane. This en- 
sures that we are indeed in the wave zone of the true 
waves where the Newtonian quadrupole formula applies. 
We have checked that /i^_ PN ex sin 2 9/r there, so that the 
gravitational wave amplitude A^q is approximately con- 
stant independent of the radius r or angle (except near 
the rotation axis where ft.+ PN vanishes). Due to the small- 
ness of /i^™ for r ^> \/(2ir), some numerical error appears 
in the cancellations of the different terms in Eq. I|35|) • That 
yields an offset in the gravitational wave signal that can 
be corrected as follows: 



2PN corrected 



= h 



2PN 



(54) 



Although the term 7 u 7i 2PN should be zero in principle, it 
is numerically comparable to /i 2PN in the wave zone; the 
parameter a corrects the offset in the waveforms. 

Notice that the link between hff N and h { ^ n is trivial in 
the near zone, where the post-Newtonian approximation 
is valid: w hf™ . We shall denote their common value 
as hj? henceforth. 



4. Initial models 

As initial models to describe either rotating neutron stars 
or rotating stellar cores at the onset of gravitational col- 
lapse we construct uniformly or differentially rotating rel- 
ativistic polytropes in equilibrium. These are obtained us- 
ing Hachi su's self-consistent field ( HSCF) method as de- 
scribed in lKomatsu et all l)l989alb|) , which solves the gen- 
eral relativistic hydrostatic equations for self-gravitating 



P = Kp 1 



(55) 



where K is the polytropic constant and 7 is the adiabatic 
index. The gauge used in the HSCF method is maximal 
slicing with quasi- isotropic coordinates (MSQI). The most 
general metric to describe these objects in the MSQI gauge 



ds 2 = ~e 2 *dP + e 2a (df 2 + r 2 d9 2 ) 



(56) 



where {i, r, 9, (p) are the coordinates in the MSQI gauge, 
and v, a, /?, and Q are metric potentials. Hereafter quanti- 
ties with a tilde are in the MSQI gauge, and all other quan- 
tities are in the ADM gauge. The hydrostatic equilibrium 
equations can be numerically integrated by prescribing a 
value for the central density p c , the rotation rate (selected 
by setting a ratio of polar radius r p to equatorial radius 
r c ), and choosing a rotation law. As it is commonly done 
in the literature we use 

as rotation law, where Q is the angular velocity of the 
fluid as measured from infinity, Q c its value at the cen- 
ter, and d = rsin0 the distance to the rotation axis. The 
positive constant A parametrizes the rate of differential 
rotation, with A — > 00 for a rigid rotator and A < r e for 
differentially rotating stars. 

For the study of rotating neutron stars, we choose 
the polytropic EOS JHEJ with K = 1.456 x 10 5 (in cgs 
units) and 7 = 2. We have constructed a sequence of 
uniformly rotating models with a central density p c — 
7.95 x 10 14 g cm -3 , and a ratio r p /r c of polar to equato- 
rial coordinate radius ranging from 1.00 (spherical model; 
labeled RNSO) to 0.65 (rapidly rotating model near the 
mass-shedding limit; labeled RNS5). A summary of impor- 
tant quantities specifying these models is given in Tabled 
Note also that models with this choice o f parameters are o f 
wi despread use in the l i teratu re (see e.g. lFont~ "all l|2002h 
or iDimmelmeier et"al1 (|2004h and references therein) and 
will be used in this work for comparison with previous 
results obtained with independent codes. 

To model a stellar iron core as initial model for simu- 
lating rotational stellar core collapse to a neutron star we 
again utilize the HSCF method with the EOS parameters 
K = 4.897 x 10 14 (in cgs units) and 7 = 4/3, chosen to ap- 
proximate the pressure profile in a degenerate electron gas. 
The initial central density is set to p c i n i = 10 10 g cm" 3 . 
Again each initial model is further determined by its ro- 
tation profile parameter A and the rotation rate, which is 
specified by the axis ratio r p /r c or equivalently by the ra- 
tio of rotational energy and gravitational binding energy, 
(i = T/\W\. . | i 

Following IDimmelmeier et~ I ll2002al) we use a hy- 
brid EOS in the core collapse simulations. This EOS con- 
sists of a polytropic part P p = Kp 1 plus a thermal part 
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Table 1. Equatorial radius r c , axis ratio r p /r c , angu- 
lar velocity 0, and ADM mass -Madm for a sequence 
of uniformly rotating neutron stars used in this paper. 
Along the sequence the rotation rate increases, from the 
spherical model RNSO to the most rapidly rotating model 
RNS5, which rotates near the mass shedding limit at 
fl K = 5.363 kHz. 



Model 


r c [km] 


r p /r c 


n/n K 


Madm [M ] 


RNSO 


14.1 


1.00 


0.00 


1.40 


RNS1 


16.2 


0.95 


0.42 


1.44 


RNS2 


17.3 


0.85 


0.70 


1.51 


RNS3 


18.7 


0.75 


0.87 


1.59 


RNS4 


19.6 


0.70 


0.93 


1.63 


RNS5 


20.7 


0.65 


0.98 


1.67 



Table 2. Parameters used in the core collapse simulations. 
The initial models are differentially rotating stellar cores 
specified by the parameter A controling the degree of dif- 
ferential rotation (cf. Eq. (JS3) and the ratio of rotational 
to potential energy, f3 = T/|W|. The collapse is initiated 
by reducing 7 from its initial value 4/3 to 71 . Additionally, 
the values for the initial equatorial radius r e and the initial 
ADM mass Madm are given. The label AxByGz of each 
model is a combination of the initial rotation parameters 
A and j3 (AxBy) and the value of 71 during collapse (Gz). 



Model 


A 







Madm 


71 




[10 3 km] 


[%] 


[km] 


[Mq] 




A1B3G3 


50.0 


0.89 


2247 


1.46 


1.31 


A1B3G5 


50.0 


0.89 


2247 


1.46 


1.28 


A2B4G5 


1.0 


1.81 


1739 


1.50 


1.28 


A4B5G5 


0.5 


4.03 


1249 


1.61 


1.28 



Pth = (7th - l)eth, where 7 th = 1-5 and e th = e - e p . 
The thermal contribution is chosen to take into account 
the rise of thermal energy d ue to shock heating . Detail s 
of this EOS can be found in iDimmelmeier et alJ l(2002al) . 
Gravitational collapse is initiated by slightly decreasing 7 
from its initial value to 71 < 4/3. As p reaches nuclear 
density p nuc = 2.0 x 10 14 g cm" 3 , we raise 7 to 72 = 2.5 
and adjust K accordingly to maintain monotonicity of the 
pressure. Due to this stiffening of the EOS the core un- 
dergoes a bounce. In Table [3 we present the main proper- 
ties of models A1B3G3, A1B3G5, A2B4G5, and A4B5G5 
used in this work to study core collapse. These models are 
identical to those with the same label s in the comprehen- 
sive cor e collapse study performed bv IDimmelmeier et alJ 
I2002b|). 

We note that when evolving the initial models con- 
structed on the basis of the HSCF method (which uses 
the MSQI gauge) with the CFC+ evolution code (which 
uses the ADM gauge), we have to consider that in general 
the two gauges differ. This could potentially lead to an un- 
suitable matching of the data describing the initial models 
when evolved with the numerical code. Let us consider the 
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Fig. 1. Deviation from conformal flatness along the equa- 
torial plane for a typical rotating stellar core initial model 
(model A1B3; solid line), and for a typical rotating neu- 
tron star in equilibrium with axis ratio r p /r e = 0.7 (model 
RNS4; dashed line). The vertical dotted line indicates the 
equatorial stellar radius for the neutron star. 



most general metric in a generic dynamic scenario, 

lij = <\>%, + hj? (ADM gauge), (58) 

% = 4>%j + 7«j> (MSQI gauge), (59) 



In general 7<y> is not transverse, so that the ADM gauge 
and the MSQI gauge differ. To quantify the differences 
between both gauges it is relevant to compare the traceless 
part of the 3-metric 7<y> and the trace 7 i - 7 7ij = 3</) 4 
itself. For an equilibrium stellar model constructed within 
the MSQI gauge one obtains 



l<ij> = (7r 



(60) 



In Fig. ^ we plot equatorial profiles of \ j rr — 7 w |/(7 y 7ij) 
for (i) the equilibrium model A1B3 used as initial data 
for core collapse simulations (lower curve; cf. Table 
and for (ii) the equilibrium rotating neutron star model 
RNS4 (upper curve; cf. Table QJ. In both cases, especially 
in the (only weakly relativistic) collapse initial model, we 
observe that deviations from conformal flatness are negli- 
gible. This fact makes the initial models built with HSCF 
method suitable for time evolution in the CFC approxi- 
mation. It also shows that the differences between both 
gauges are very small, namely of the same order of magni- 
tude as 7<y> ■ Consequently, the use of initial models com- 
puted in the MSQI gauge for evolutions using the ADM 
gauge is entirely justified since the differences are negligi- 
ble and only appear at the second post-Newtonian order. 

5. Rotating neutron stars 

5.1. Effects of the boundary on the metric solution 

As mentioned in Sect. 13.21 the location where the bound- 
ary conditions for the CFC+ intermediate potentials are 
imposed can have a significant influence on the accuracy 
of the metric solution. Additionally, the extension of the 
numerical grid is also of paramount importance for the 
numerical computation of the CFC equations for the con- 
formal factor, the shift vector, and the lapse function. As 
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r [km] 

Fig. 2. Effects of the outer boundary conditions on the 
equatorial radial profile of the intermediate potential S r 
(top panel) and the metric component hJ^F (bottom panel) 
for the rotating neutron star model RNS4. Zero-value 
boundary conditions are imposed at 62 km (dashed line), 
at 10 4 km (dashed-dotted line), at 10 6 km (dashed-dot- 
dotted line), and at 10 11 km (solid line). Alternatively, 
monopolc boundary conditions are imposed at 62 km 
(open boxes) and at 10 6 km (filled circles). Overplotted 
to S r is the monopole behavior (dotted line). 



the sources of these equations have non-compact support, 
it is necessary to integrate out far enough in radius to 
obtain the desired accuracy of the numerical solution. 

To accomplish this we use an extended (vacuum) grid 
going far beyond the actual stellar boundary. Monopole 
behavior at the outer boundary of this vacuum domain 
has been checked by comparing with calculations done 
imposing zero values for the potentials as outer bound- 
ary conditions at extremely large distances (see Fig. [SJ. 
Convergence tests with different parameters for the outer 
grid (number of cells, distance of the outer boundary, am- 
plification factor of the grid) have been performed. The 
conclusion of these tests is that for simulations of rotating 
neutron stars we need to use a numerical grid consisting 
of ne = 30 angular cells and n r — 250 radial cells (of 
which 100 are equally spaced inside the neutron star and 
150 are logarithmically spaced for the atmosphere) to cor- 
rectly capture the conformally flat part of the metric. In 
order to impose outer boundary conditions for the CFC+ 
metric potentials, an extra grid of 70 radial cells extending 
out to 10 6 km needs to be added. 

5.2. CFC+ corrections to the 3-metric 

We turn now to measure the magnitude of the CFC+ cor- 
rections to the metric of our sample of rotating neutron 
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Fig. 3. Radial profiles at the equator of the non- vanishing 
components of hj? for the sequence of models RNS0 to 
RNS5. The strength of the correction hj? increases with 
the rotation rate. The equatorial radius of each model is 
given in Tabled 



star models (cf. Table First, from the distribution of 
the hydrodynamic variables in the equilibrium models pro- 
vided by the HSCF method we recompute the CFC+ met- 
ric components. For equilibrium configurations, i.e. in an 
axisymmetric stationary spacetime, the components hJJ 
and hJJ vanish. All other components are shown in Fig. [3] 
for different models with increasing rotation rate. Note 
that as discussed in Sect. 0] the ADM gauge is not quasi- 
isotropic in the CFC+ approach, because hJJ ^ and 
^ Kqq . Therefore, a direct comparison with the full 
general relativistic metric calculated in the MSQI gauge is 
not possible. As expected, hjj vanishes for the spherical 
case RNS0, in which both CFC and CFC+ agree exactly 
with general relativity, and there is no traceless and trans- 
verse part. As we increase the rotation rate from RNS1 to 
RNS5, the value of hj? increases, resulting in a depar- 
ture of the metric from conformal flatness. As KF7' is a 
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Table 3. Frequencies of small-amplitude quasi-radial pul- 
sations for model RNS4. We compare the frequencies ob- 
tained from simulations with the present code (using ei- 
ther the CFC or the CFC+ approximation) with those 
obtained independently from a 3D full general relativistic 
code (GR). The results are extracted from time evolutions 
where the spacetime metric is kept fixed (Cowling approx- 
imation). The relative differences between the CFC+ and 
the GR code are shown in the last column. 



Mode 


jCFC 


jCFC+ 




Rel. diff. 




[kHz] 


[kHz] 


[kHz] 


[%] 


F 


2.48 


2.48 


2.468 


0.5 


HI 


4.39 


4.39 


4.344 


1.1 


H2 


6.30 


6.30 


6.250 


0.8 



correction to the conformally flat 3-metric (which is of or- 
der unity, 4> ~ 1), it can be shown that even in our most 
extreme case RNS5, which is very near to the mass shed- 
ding limit, the values of hj^ amount to a correction of 
only about 1%. Therefore, qualitative differences with re- 
spect to the CFC approximation are not expected in the 
dynamics of such objects. 

5.3. Oscillations of rotating neutron stars 

A further test of the new approximation for the metric is 
provided by the study of pulsations of rotating neutron 
stars. To this aim we perturb the neutron star models de- 
scribed in Sect. 0] with a density perturbation of the form 
Ppert = P [1 + acos(7rr/(2r c ))], where a is an arbitrary 
parameter controlling the strength of the perturbation. 
The perturbed models are evolved in time in two different 
ways, either considering coupled evolutions for the hydro- 
dynamics and the metric, or evolving only the hydrody- 
namics in a fixed background metric corresponding to the 
metric provided by the elliptic solvers at the first time step 
(an approximation commonly referred to in the literature 
as the Cowling approximation). Both metric approxima- 
tions, CFC and CFC+, are used to compare the respective 
results. 

The oscillations of the stars can be observed in vari- 
ous hydrodynamic and metric quantities. In particular we 
monitor the central density p c , the radial velocity v r at 
half the stellar radius, and the gravitational wave ampli- 
tudes A^o extracted with the quadrupole formula. When 
Fourier transforming the time evolution of these quanti- 
ties, distinctive peaks appear at the same (discrete) fre- 
quencies for any of these variables. Those frequencies can 
be identified with the qua si-normal modes of pulsation 
of the star, as described in iFont et n ,lJ feOOfjl . To further 
validate our approach we present the quasi-normal modes 
calculated in the CFC and CFC + approximation in com- 
parison with those reported by IFont et ahl lj2002j) , which 
are calculated with a 3D code in full general relativity 
(without any approximation). 



Table 4. Fundamental mode frequency /p of small- 
amplitude quasi-radial pulsations for a sequence of rotat- 
ing polytropes with increasing ratio of polar to equatorial 
radius r p /r c . We compare the frequencies obtained from 
simulations with the present code using either the CFC or 
the CFC+ approximation with those obtained indepen- 
dently from a 3D full general relativistic code (GR). The 
results are extracted from coupled spacetime metric and 
hydrodynamic time evolutions. The relative differences be- 
tween the CFC+ and the GR code are shown in the last 
column. 



Model 


r p/ r c 


i-CFC 
JF 


J-CFC+ 
Jf 

[kHz] 


f GR 
JF 


Rel. diff. 






[kHz] 


[kHz] 


[%] 


RNS0 


1.00 


1.43 


1.43 


1.450 


1.4 


RNS1 


0.95 


1.40 


1.40 


1.411 


0.8 


RNS2 


0.85 


1.34 


1.34 


1.350 


0.7 


RNS3 


0.75 


1.27 


1.27 


1.265 


0.4 


RNS4 


0.70 


1.24 


1.24 


1.247 


0.6 


RNS5 


0.65 


1.21 


1.21 


1.195 


1.0 



Table [3] shows the mode-frequencies for the fundamen- 
tal mode (F), as well as for the first (HI) and second (H2) 
harmonics obtained in evolutions of model RNS4 in which 
the spacetime metric is kept fixed (Cowling approxima- 
tion). The accuracy in the frequency values depends on 
the total time of the evolution, increasing as the evolution 
is extended. We have evolved all models for 30 ms, finding 
no significant deviations in the hydrodynamic profiles with 
respect to the original profiles of the equilibrium models. 
For such an evolution time the FFT yields a maximum 
frequency resolution of 0.03 kHz. Table shows that no 
differences can be observed between the mode-frequencies 
computed with CFC and CFC+, and that, in addition, 
there is very good agreement with the general relativistic 
results since the reported values are within the affordable 
resolution in frequency. 

The corresponding results for the case of coupled evo- 
lutions of the spacetime metric and the hydrodynamics are 
shown in Tables 01 and for the fundamental mode and 
the first harmonic, respectively. In these simulations, for 
the sake of computational efficiency and without affecting 
the dynamics, the metric is calculated every 100th hydro- 
dynamic t ime steps and extrapolated in between, as ex- 
plained in iDimmelmeier et alJ l|2002al) . As in the Cowling 
simulations, all models are evolved for 30 ms. Even for 
the more rapidly rotating models, no differences in the 
frequencies from the CFC and CFC+ simulations can be 
found. This result can again be explained by the small- 
ness of the components of hj^- , which does not modify 
the dynamics considerably. Furthermor e, the result s agree 
to high precision with the GR results of lFont et al.l l|2002t) 
within the limits set by the temporal and spatial resolu- 
tion. 

We emphasize that, for accurately extracting the os- 
cillation mode frequencies, the code has to maintain the 
initial equilibrium configuration in a hydrodynamical evo- 
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Fig. 4. Time evolution of hydrodynamic and metric quantities for the regular collapse model A1B3G3 (left) and the 
rapid collapse model A1B3G5 (right). The top panels show the central density p c (solid line) and lapse function a c 
(dashed line). Both the CFC and the CFC+ results overlap. Nuclear matter density p nuc is indicated by the horizontal 
dotted line. The second panels from the top display the relative difference a of p c (solid line) and a c (dashed line) 
between the simulation using CFC and CFC+. In the third panels the CFC+ evolution of the absolute value of hJJ 
(solid line), hJJ (dashed line), /iJJ (dashed-dotted line), as well as the trace of hj? (dotted line) is shown, all measured 
at the center of the star. Note that hJJ and hTE cannot be discerned, as they practically overlap. The bottom panels 



depict the evolution of the maximum absolute value of h^g (solid line), (dashed line), and hg^ (dashed-dotted 
line), respectively. Again the latter two quantities closely coincide. The vertical dotted line in all panels marks the 
time of bounce %■ 



Table 5. Same as Table 0] but for the frequency /hi of periods). Independently of the approximation assumed for 

the first harmonic mode. In the model RNS5 this harmonic the metric (either CFC or CFC+) and the (small) gauge 

was not sufficiently excited by the perturbation chosen for mismatch, we have tested that our code is able to perform 

a clear identification of its frequency. that task successfully. 



Model 


r p /r c 


/•CFC 

/hi 


J-CFC+ 
JH1 


rGR 

/hi 


Rel. diff. 






[kHz] 


[kHz] 


[kHz] 


[%] 


RNS0 


1.00 


3.97 


3.97 


3.958 


0.3 


RNS1 


0.95 


3.87 


3.87 


3.852 


0.5 


RNS2 


0.85 


3.95 


3.95 


3.867 


2.0 


RNS3 


0.75 


3.98 


3.98 


4.031 


1.3 


RNS4 


0.70 


4.02 


4.02 


3.887 


2.0 


RNS5 


0.65 






3.717 





lution for many rotation periods (usually several tens of 



6. Rotational core collapse 

We now present results of simulations of rotational core 
collapse to neutron stars. The core collapse models we 
have selected (see Table |2J| are representative for the dif- 
ferent types of collapse dynamics and gravitational ra- 
diation waveforms observe d in the CFC simulations of 
iDimmelmeier et"al1 l|2002a|) : A1B3G5 as type I (regular 
collapse), A2B4G5 as type II (multiple bounce collapse), 
A1B3G5 as type III (rapid collapse), and A4B5G5 as a 
case with extreme rotation, i.e. a strongly and highly dif- 
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Fig. 5. Same as Fig. 0] for the multiple bounce collapse model A2B4G1 (left) and the extremely rotating collapse 
model A4B5G5 (right). Both models centrifugally bounce before reaching p nuc . Note that for model A4B5G5 we plot 
the maximum density /9 max instead of p c , as this model has a toroidal density configuration. 



ferentially rotating core with an initial torus-like structure 
which is strongly enhanced during collapse. We use a nu- 
merical grid consisting of ng — 30 equally spaced angu- 
lar cells and n r = 300 logarithmically spaced radial cells 
covering the star. In order to calculate hj^- and extract 
waveforms, an extra grid of 300 radial cells extendig out to 
10 11 km needs to be added (600 cells for model A2B4G1). 

6.1. Collapse dynamics 

In Figs. 0] and 03 we compare the time evolution of se- 
lected matter and metric quantities for all four collapse 
models considered using both the CFC and the CFC+ 
approximation. We show the time evolution of the cen- 
tral density p Cl which is a representative quantity of the 
hydrodynamic evolution. We present p c for all models ex- 
cept model A4B5G5 for which the time evolution of the 
maximum density p max is used instead. In this model the 
density maximum is not attained at the center due to 
the strong differential rotation. The evolution of the cen- 
tral density in models A1B3G3 and A1B3G5 (see Fig. EJ) 
shows a distinctive rise during collapse until p c reaches 
its maximum at the time of bounce tb (at tb ~ 49 ms 



and tb ~ 30 ms, respectively). Later on the density os- 
cillates around the new equilibrium value of the compact 
remnant (which can be identified with the new-born proto- 
neutron star). These oscillations are highly damped due to 
the presence of an extended stellar envelope surrounding 
the proto-neutron star. Note that in models A2B4G1 and 
A4B5G5 (see Fig. 01 the collapse is stopped by rotation be- 
fore nuclear matter density is reached, as strong centrifu- 
gal forces build up during the collapse. Hence, the evolu- 
tion of model A2B4G1 is characterized by the presence of 
consecutive multiple bounces, while the centrifugal hang- 
up in model A4B5G5 causes a single bounce below nuclear 
matter density, leaving a low density proto-neutron star 
behind. 

The top panels of Figs. 01 and El also show the central 
lapse function a c (dashed line; labels on the right vertical 
axis) . In the CFC+ approach the new hj? terms directly 
couple to the metric equation 129fl for a, while they couple 
indirectly to the metric equations for <jj and [3 l through a 
itself. The evolution of the lapse closely mirrors that of the 
density, decreasing while the density increases, i.e. while 
the star contracts, and vice- versa. 
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For the four collapse models considered here, there is 
no direct visual evidence of discrepancies between the CFC 
and CFC+ results. The corresponding curves for CFC and 
CFC+ in the top panels of Figs. 0] and coincide per- 
fectly in the case of both the lapse function and the den- 
sity. Therefore, as no appreciable differences are visible, 
we plot in the second panels (from the top) of these fig- 
ures the relative differences a of these two quantities be- 
tween CFC and CFC+. Maximum differences of the order 
of ~ 1% are found in the density evolution (solid line) for 
the strongly differentially rotating models A2B4G1 and 
A4B5G5 (see Fig. In the two other models the dif- 
ferences are two orders of magnitude smaller. The lapse 
function (dashed line) shows even smaller differences be- 
tween CFC and CFC+, the maximum values of a being 
smaller than 0.1% even for the rapidly rotating models. 

The time evolution of the diagonal components of 
hj^- at the center (r = 0) are also plotted in Figs. ^ 
and (third panels from the top) , along with the trace 
of hJ 3 T . Correspondingly, the maximum absolute values of 
the off-diagonal terms of hj^ are displayed in the pan- 
els at the bottom. As expected, the various components 
of h J T arise and increase when deviations from sphericity 
occur. The profiles show that in all collapse simulations 
hjj 1 - is quite small in comparison with the isotropic part 
which is of order unity. It can be seen that models with 
strong gravity but small asphericities (as A1B3G3) and 
models with weaker gravity but more apparent deviations 
from sphericity (as A2B4G1 or A4B5G5) all reach values 
for hJ } T of similar magnitude. Note that the components 
hJJ and hj£ rapidly decrease after the bounce, because 
a quasi-equilibrium configuration is reached in the new- 
born proto-neutron star. In all cases considered the trace 
of hj^~ is much smaller than the hj? components them- 
selves, i.e. numerically hj^ is traceless to high accuracy, 
and also remains traceless during the entire evolution. In 
addition, we have checked the transverse character of hj^ , 
i.e. V l /i^ T = 0. The latter expression is found to be com- 
patible with zero, as the dimensionless quantity rV'/iJ T 
is much smaller than hj^ . 

The radial profiles of hj^ are very similar for all col- 
lapse models we have analyzed except for model A4B5G5 
that collapses off-center, with a torus-like structure. In 
Fig. we compare this model with a model in which 
the maximum density is reached at the center (A1B3G3). 
The profiles are depicted at the instant of maximum den- 
sity (t b = 48.9 ms for model A1B3G3 and t b = 31.2 ms 
for model A4B5G5) and at the end of the simulation, 
when the system has reached an equilibrium state. For 
the spheroidal model A1B3G3 the maximum values of 
hj^ are reached at the center, and the components hJJ 
and /iJJ have local maxima inside the star. However, in 
the toroidal model A4B5G5 the maximum values are off- 
centered, while the three components exhibit their peak 
value inside the torus. Note that the strong deviations 
from sphericity generate in model A4B5G5 larger values 
of hj^ than in model A1B3G3 at the time of bounce, but 



once the torus collapses to the final oblate star, the values 
become smaller than for the regular collapse model. 

Table summarizes the results of all collapse simula- 
tions, including relevant information to calculate the size 
of the near zone \/(2ir) needed for the gravitational wave 
extraction which we discuss next. 

6.2. Gravitational radiation waveforms 

Gravitational waves from the collapse simulations dis- 
cussed in the preceding section have been calculated for 
both the CFC and the CFC+ approximation of the field 
equations, using the quadrupole formula. In addition, in 
the case of CFC+ they have also been extracted directly 
from hj? evaluated in the wave zone. The radial exten- 
sion of the near zone X/(2tt) can be calculated from the 
approximate size of the source r e and the typical timescale 
of gravitational wave emission l/f max . Results for each 
collapse model are listed in Table El which also gives the 
distance r cx at which the waveforms are actually extracted 
from hJ 3 T . For all models this distance is much larger than 
X/(2n) (by a factor of about 50), which ensures that the 
wave extraction is done in the wave zone far away from 
the sources. 

The gravitational waveforms are displayed in the top 
panels of Figs. and |H1 respectively. The waveforms ex- 
tracted using the quadrupole formula (solid lines) are very 
similar for CFC and CFC+ and would not be discernible 
in the figures. Thus, only the CFC+ waveforms are plot- 
ted, along with the absolute differences cr a b s with respect 
to CFC (shown in the bottom panels of each figure). For 
all collapse models considered the differences are smaller 
than 0.1% of the signal maximum. This is expected from 
the fact that the quadrupole formula involves an integral 
of hydrodynamic quantities, and from the observation that 
the modifications in the collapse dynamics between CFC 
and CFC+ are not significant, as mentioned in Sect. 16. ll 

Concerning the waveforms extracted directly from hj^- 
(dashed and dotted lines in Figs. and it can be 
seen that when directly using Eq. (|52|l to calculate the 
wave amplitude A^q , the resulting signals are larger at 
bounce for all models. After bounce these signals show 
an offset for the models with stronger gravity (A1B3G3 
and A1B3G5, see Fig.UJ), where the waveform amplitude 
should actually approach zero, because the pulsations of 
the new-born proto-neutron star are rapidly damped and 
the star tends towards an equilibrium state. If the signals 
are corrected by means of Eq. I|54[l , the offset disappears 
and the gravitational wave amplitude agrees remarkably 
well with the waveforms calculated with the quadrupole 
amplitude. Although the extraction methods are not re- 
ally independent, the agreement found between the two 
ways of calculating waveforms in the CFC+ approach is a 
consistency check for the calculation of the hj^ , because 
the asymptotic behavior given by Eq. (|52J) can be assessed 
numerically this way. 
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Table 6. Summary of various quantities that characterize the different core collapse models. The table shows the time 
of bounce t\>, the maximum density at bounce pb, the maximum density reached after the bounce pf , the gravitational 
wave amplitude at bounce |^4fp | b as measured using the quadrupole formula, the dominant frequency of oscillations 
of the proto-neutron stars / maX) the radius of the star after bounce and ringdown r c f (defined as the radial location 
along the equator where the density first falls below 10% of the maximum density), the size of the near zone X/(2tt), 
and the distance r cx at which gravitational waves are extracted from hJ T . 



Model 


t b 


Pb 


Pf 


1 4 E2 I 

1^20 | b 


/max 


7ef 


A/(2tt) 






[ms] 


[10 14 g cm" 3 ] 


[10 14 g cm" 3 ] 


[cm] 


[Hz] 


[km] 


[km] 


[10 3 km] 
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Fig. 6. Radial profiles of hj? at the equator for model A1B3G3 (left) and model A4B5G5 (right) at the time of 
maximum density (upper panels) and at the final time of the simulation (lower panels). The curves plotted correspond 
to hJJ (solid line), KgJ (dashed line), (dotted line), and the trace of tif? (dashed-dotted line), respectively. 



The agreement holds for all cases except for 
model A2B4G1, where the amplitudes obtained by the 
quadrupole formula and from the hj^ differ by about 50% 
both at the main bounce and at the subsequent bounces 
even after the offset correction is applied. Note that for 
this particular model the size of the near zone is very 
large, A/(27r) = 884 km. Therefore, an accurate extrac- 
tion of the waveforms can only be performed at a radius 
very far away from the star. We thus set the extraction 
radius r cx = 4 x 10 4 km. As a consequence the extended 
grid needs to be covered with at least 600 radial zones in 
order to avoid too extreme logarithmic cell spacing, which 
would be the source of numerical inaccuracies in the fall- 
off behavior when solving the elliptic equations with finite 
difference methods. 



7. Conclusions 



We have presented a new approximation for the Einstein 
field equations, which we call CFC+. We have tested its 
suitability for simulations of rotating neutron star space- 
times, both for equilibrium models and for configurations 
which are formed after gravitational core collapse. This 
approach is based on second order post-Newtonian cor- 
rections from conformal flatness, i.e. CFC+ represents an 
extension of the CFC (or Isenberg- Wilson-Mathews) ap- 
proximation. The derivation of the extended system of 
equations has been presented in great detail, as well as 
the boundary conditions to apply when numerically solv- 
ing them. All CFC+ field equations have elliptic charac- 
ter, since at second post-Newtonian order the hyperbolic 
character of the Einstein equations disappears. This is a 
consequence of the fact that the time derivatives of 
appear first at the 2.5th post-Newtonian order. 
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Fig. 8. Same as Fig.Hfor the multiple bounce model A2B4G1 (left) and the extremely rotating collapse model A4B5G5 
(right). 



We note in passing that solving elliptic equations en- 
sures the numerical stability of the solution and avoids the 
numerical problems sometimes encountered in long-term 
evolutions of strongly gravitating systems when using the 
3 + 1 formulation of general relativity. On the other hand, 
the price to pay for using this approximation is that grav- 
itational radiation reaction on the dynamics, caused by 



gravitational waves carrying away energy and angular mo- 
mentum from the system, is absent. However, in the case 
of models where a comparison of our results to fully gen- 
eral relativistic results is possible, we have checked that 
the absence of gravitational back-reaction does not affect 
the results significantly. In scenarios such as the merging of 
compact binaries (not investigated here), this effect would 
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indeed be important, but only at late times. Hence, CFC+ 
should also be a good approximation to model phenomena 
occurring on dynamical timescales such as the final stages 
of the plunge and merger. 

We have compared the n ew CFC+ approxim a tion with 
the CFC approach used bv iDimmelmeier et alJ l|2002alb|) 
in two different scenarios, oscillating relativistic stars and 
core collapse to a neutron star. In the case of pulsating 
neutron stars, we find that there are no differences in the 
calculation of the quasi-radial normal mode frequencies of 
those objects, even in the most extreme situations consid- 
ered when the star is rotating at the maximum allowed 
rate (i.e. at the mass-shedding limit). It has been pos- 
sible to compare our results directly with fully general 
relativistic computations and, again, a very close agree- 
ment is found. Furthermore, our simulations of stellar 
core collapse to neutron stars covered well the basic mor- 
phology and dynamics of c ore collapse types studied by 
IDimmelmeier et"al1 l)2002tJl . including the extreme case 
of a core with strong differential rotation and torus-like 
structure. Once more, no significant differences between 
the two approximations are observed. Therefore, we can 
conclude that second post-Newtonian corrections to CFC 
do not significantly improve the results when simulating 
the dynamics of core collapse to a neutron star as well as 
when investigating the evolution of neutron stars in equi- 
librium. 

Regarding the gravitational wave extraction we have 
not observed any substantial differences between CFC and 
CFC+ as well. The comparison has been carried out using 
the quadrupole formula, commonly employed in the liter- 
ature to extract gravitational waveforms. In addition we 
have also calculated the gravitational waves directly from 
the hj? , which permits a straightforward use of the space- 
time metric to study the gravitational wave generation 
mechanism from the near zone to the wave zone. Although 
the waveforms computed with the latter approach cannot 
be regarded as an independent way of calculating gravita- 
tional wave signals, it nevertheless provides a good consis- 
tency check of the CFC+ approximation that has served 
to validate the numerical scheme we use to calculate h? T . 

The main conclusion of this work is the assessment 
of the CFC approximation as a highly suitable alterna- 
tive to the full Einstein equations in axisymmetric scenar- 
ios, involving rotating neutron stars in equilibrium and 
as end states of core collapse. These findings are sup- 
ported by two facts: First, we have demonstated that sec- 
ond post-Newtonian corrections do not play an important 
role in neither the dynamics nor the gravitational radia- 
tion waveforms of core collapse. This suggests that higher 
order post-Newtonian corrections will be completely in- 
significant at least on dynamic timescales. Secondly, the 
direct comparison of the CFC approach with exact fully 
general relativistic simulations of pulsating neutron stars 
yields normal-mode frequencies in excellent agreement. 
Furthermore, comparisons of the CFC approach with fully 
general rel ativistic simulations have als o recently been re- 
ported bv lShibata &; Sekiguchil l|2004a|) in the context of 



axisymmetric core collapse simulations. Again, the dif- 
ferences found in both the dynamics and the waveforms 
are minute, which highlights the suitability of CFC (and 
CFC+) for performing accurate simulations of those sce- 
narios without the need for solving the full system of the 
Einstein equations. 

With this investigation we have shown that a numeri- 
cal code based on CFC is a very useful tool to investigate 
core collapse to neutron stars in general relativity. This 
approach is suitable to form the basis of a future core col- 
lapse supernova code which can be gradually extended in 
various directions to incorporate additional physics such 
as realistic equations of state, magnetic fields, and eventu- 
ally neutrino transport. In the near future we plan to fur- 
ther validate the CFC+ equations in other scenarios where 
higher densities are present (e.g. collapse to a black hole), 
as well as in the fully three-dimensional case (namely to 
investigate the merging of neutron stars). Such scenarios 
are in principle beyond the range of applicability of the 
CFC approximation, but can possibly still be handled in 
a satisfactory way with the new CFC+ approach presented 
in this work. 
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Appendix A: Metric equation for the lapse 
function in CFC+ 

The correction to the CFC metric up to second post- 
Newtonian order may be obtained in the Hamiltonian 
framework of ADM within the eponym gauge (see also 
iRegge fc Teitelboiml ill 9741) for an ana lysis in asymptot- 
ically flat spacetimes and Irlolml £.985) for a generaliza- 
tion to isolated perfect fluids). The original canonical vari- 
ables are chosen to be the 3-metric jij and its conjugate 
momentum cV j / (16ttG) with ir ij = -^/j(K ij ~ K-y ij ), 
but once the coordinate system has been fixed, four of 
the six remaining field degrees of freedom are eliminated 
by imposing the Hamiltonian constraint l|14|) and the 
three momentum constraints <|15|) . These four field de- 
grees of freedom correspond to the conformal factor <f> 
and to the symmetric trace-free part tt^l °f the tensor 
2A- 1 (7 lfc V fc V i 7r^) - A- 2 (f"7^V fc ViV m V„7r fci )/2, re- 
spectively. 

Only two transverse trace-free (TT) field variables are 
left, namely hj^ — Jij — y'Jij on the one hand, and 7r^ T = 
tt 11 — 7Tll on the other hand. By construction, we have: 

7 fci V fe ^ T =0, f^ T = 0, (A.l) 

Vj71Vj? T = 0, 7ij 7r TT ~ 0. (A. 2) 
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The reduced Hamiltonian H is obtained by substi- 
tuting <j> with cj)[hpj , ir^? T ] and with 7r lJ '[h-J, ir^] 
in the Hamiltonian of general relativity for asymptot- 
ically flat spacetimes. The contributions of the super- 
Hamiltonian and super-momentum densities vanish, so 



that H [matter variables, h. 



XT 



7t T t] is given by the surface 



integral defining the ADM mass "on a shell" , 



H = 



WttG 

„4 



2ttG 



d 2 S l Jf? k (VjTifc 
d 3 a; y^Acj), 



(A.3) 



as rcf. lSchwingerllll963ln 

fR = - 7 2 V 4 V, (ry y ) 



( 77 ^) V,-7Vj7 



7 (77J 7 ) Vj (77^) Va k i 



+ 2 (^ lfc ) (^) V^ViTiM 
+ ^7 (77*) V fc (77*) V i7j 7 



4 177 



77 lfc ) Vfc7V i 7 i ;. 



(A.7) 



where the surface element dS l V7 refers to the flat metric. 
The reduced Hamiltonian (|A.3(1 contains the full dynami- 
cal hn^mMioji_aboutthe_sj^ em. In particular, as shown 
bv lReeee fc Teitelboiml l|l974|h the field evolution is gov- 
erned by the equations of motion 



with 



-TTfcZ 
hj 



16ttG TTfc; SH 

167rG .jj Si? 
7a 7TTfc; 



(A.4) 
(A.5) 



(A.6) 



and a similar formula for T^Tfcz- The role of these op- 
erators is to ensure the transverse trace-free projection 
of the Frechet derivative 8H/8ir^ T and SH/ShJ^- , respec- 
tively. The calculation of H in terms of D* , S* , P, as well 
as the field variables hj^ and 7r^ T , can be done essen- 
tially by eliminating <j) in Eq. (|A.3(I with the help of the 
Hamiltonian constraint (|14|l . This is achievable in pertur- 
bative treatments such as the post-Minkowskian or the 
post-Newtonian ones, consisting of the formal expansion 
of all quantities at play in powers of the gravitational con- 
stant G or of the inverse of the square of the speed of light 
1/c 2 . 

In the course of eliminating <f>, we use Eq. I|14|) in 
a more explicit form. For this we first express the 3- 
curvature R as a function of <f> and hj^ . By making ex- 



jk — 



llij 



and 



tensive use of the relations 7^ V/7 
Vfe7 = (77 lJ )Vfc7y, the combination j 3 R can be written 



By definition, the determinant 7 it is equal to the an- 
tisymmetric sum of products 3! 7 p ' 1 7 2 -7 3 ' r '7pi7 9 27r3 = 
7 P ' t 7 J -7 fe ' r 7pi7gj7rfe on a Cartesian grid, the square brack- 
ets denoting antisymmetrization of non-underli ned in- 
dices. Its explicit expression is given, e.g., by ISchaferl 
l|l985|) . Similarly, the ij-component of the comatrix 77 1 -? 



is known to be e (£ r mn7"' r 7 mp 7 Tl<? 7pfc7 9 /)/2, where e 
the permutation operator. From the identity s l e 



ikl 



IS 



j in n 



3! S^S^Sn it is straightforward to obtain jikjjijj 1 * 1 as a 
function of <f> and hj^ . After inserting the relations for 
7 and 77 y in the right-hand side of Eq. I|A.7|) . it is ex- 
panded in powers of hj^ and truncated consistently at 
the post-Newtonian level of hj^hj^-, denoted as 0(h 2 ). 
At this point we have an expression for jR as a function 
of A(f>. 

In addition to R, there is a second contribution to 
the left-hand side of the Hamiltonian constraint surviv- 
ing in the absence of matter, namely K^K 1 ^ — K 2 = 
[•7rj-7r| — ^.{^i) 2 ]/l- In the ADM formulation, the momen- 

7r^ T . The first term twL 



turn ir lJ decomposes into ~k %3 



ll 

is of order 1/c 3 , being a sum of derivatives of Poisson in- 
verse operators acting on Vj-zr 4 - 7 — C(V jK l i), which is it- 
self C(l/c 3 ) according to the momentum constraint equa- 
tion. The second term, linear in hJ 3 T , is transverse and 
trace-free, hence 7ij7TxT = 0- Moreover, the ADM gauge 
condition implies %j~K %3 — Jijir^ = 0. 

Finally, we consider the matter source term E in the 
Hamiltonian constraint. The corresponding density E* = 
y/^E may be written as 



16ttG 



-E* 



IGirG 



4-mn ^m^n 



1/2 



2 7 7 n kl c2 ^ 2 



"PI 



(A.- 
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This yields an elliptic equation for V = 2(0—1) = 
C(l/c 2 ) up to 0(h/c e ) and 0{h 2 /c 2 ) corrections: 



-4AV = ~4:V i (j)V j (l)f k ^ jl h^ 
1 

+ 4 



I ) 




LL 7T XX -f- 7T TT 7r TT 



(D*h) 2 + ^~ 4 <y mn -lS 



1/2 



Sp 



(A.9) 



As the terms containing a factor hJ 3 T or 7r^ T are propor- 
tional to the coupling constant G/c 2 of general relativity, 
Eq. EE reduces to AV~ = -47rGL>7c 2 + C(l/c 4 ) at the 
lowest post-Newtonian approximation. Thus, if we intro- 
duce the Newtonian potential U defined as the smooth 
solution of the Poisson equation l|28(l vanishing at spa- 
tial infinity, we have V = U/c 2 + C(l/c 4 ), plus a possi- 
ble harmonic function. Assuming an asymptotically flat 
space-time, this function must tend asymptotically to- 
wards zero while being regular, and so it has to be iden- 
tically zero. Another important piece of information pro- 
vided by Eq. I|A.9|) is the value of the lowest order contri- 
bution to the potential V that depends on the field vari- 
ables. It is given by the equation 



pure matter part, 



(A.10) 



which shows incidentally that <f> = 1 + V/2 is not af- 
fected by a non-zero hj^ at the leading post-Newtonian 
approximation in the field. Inserting the resulting expres- 
sion for the conformal factor into the Hamiltonian con- 
straint 1A.9I). we arrive at 



1 / - - S*S 
AAV = — j 2ViUVjU + 8irG- 



D* 



7 7 h kl 



1 



+ 77 fc ¥ m 7 Jn V fe /^V^T 



4 



kl 
TT 



+ total derivative + pure matter part. (A. 11) 

The pure matter part has been kept aside because it does 
not enter the computation of hj^ . The total derivative 



terms are irrelevant for the investigation of the field evo- 
lution. Indeed, by virtue of relation l|A.3|) . the Hamiltonian 
is given by the spatial integral of —c A AV/ (AnG). Thus the 
dynamics of the gravitational interaction is described by 



H field = 

+ mt. 15 



D* 



1 



O 



4 

+ liklji 
h 



kl 
TT 



o 



(A.12) 



in agreement with lSchaferl l|l990|) . The Hamilton equations 
provide the evolution of the field. They take the explicit 
form 

Othij = 2c lij [ikmlln^LL + 7T TT)\ 

■o(\) + o( k - 



2ciikijin, 



kl 



O 



O 



(A.13) 



D* 



ran ) 



o 



-A ( 7 fcm f n h 



O 



(A.14) 



The non-conformally flat part of the 3-metric appears first 
at the second post-Newtonian approximation. Its leading 
order is obtained by inserting the above expression for 
9 t 7r^ T into the time derivative of Eq. I|A.13I) . The result- 
ing equation is of wave type. In the near zone, all terms of 
order 1/c 6 may be neglected, in particular the time deriva- 
tive contribution to the dAlembertian operator. Hence 
the non-conformally flat part of the 3-metric satisfies 



-TTfcZ 
hj 



AV k UViU +16nG SkSl 



D* 



O 



(A.15) 



which is identical up to second post-Newtonian order to 
Eq. J2EJ- 

The equation for the lapse function a is derived from 
the gauge conditions (|17II18|) combined with the evolution 
equation (|13|) for the extrinsic curvature K%j. From the 
identity n^iij = 0, we deduce that the trace is negligible 
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at this level of approximation: 



f J K k 



2V7 



O 



1 



(A.16) 



By contracting Eq. (|13|) with the inverse 3-metric 7 y , we 
arrive at 

^d tKij = (TK^+O (A) = O (i) , (A.17) 

from which follows that 



V;V*a 
4ttG 



+ ai? 
t (-5 + 3£0 



(-2aiQ fe + 2V t [3 k ) 

o- 1 



(A.18) 



Due to the fact that K lk is symmetric and trace-free 
modulo 0(h/c 3 ) corrections, the product K lk (—2aKi k + 
2V J( 9fc) can be written as K tk (-2aK <ik> + 2V<i/3fe>) + 
0(h/c 6 ). The terms inside the parentheses vanish accord- 
ing to the symmetric trace- free version of Eq. Ijl2(l , 

-2aK <l3> + 2V <i /3 j > = - c d tl<ij> = O (j^J , (A. 19) 

so that K lk (-2aK lk + 2V l /3 fc ) = C(l/c 8 ). Next, we see 
from the Hamiltonian constraint equation that the inter- 
action and field parts of the scalar curvature R appearing 
in Eq. |TA~18|) are actually of order 0{h/c 4 ) = C(l/c 8 ). 
On the other hand, we know that E* = 0{hjc^) + 
pure matter terms. We have similar equalities for E and 
S = Si. Therefore, we find 



pure matter part 



= Aa - ViVjaf k j jl hl?. (A.20) 



negligible at the second post-Newtonian level, and will not 
be computed here. As already pointed out, the equation 
for the conformal factor remains unaffected at that level 
as well. Thus in general, all CFC equations except the one 
for the lapse function remain unaltered. 



Appendix B: Inversion of the equation for the 
3-metric correction in CFC+ 

To express the second post-Newtonian correction hjj 1 in 
CFC+ with the help of intermediate potentials and thus 
to invert Eq. H26[) into Eq. I|35ll in Sect. 12.41 we proceed in 
three stages: 

(i) We make the action of 7^™ explicit in Eq. (|2*uT) . 
The result is integrated formally by means of the Poisson 
integral operator A -1 . By virtue of its linearity prop- 
erty, we obtain a weighted sum of Poisson potentials 
of generic type A~ 1 F mn (up to possible index contrac- 
tions), or super-potentials of the form A _2 VjVj.F mn — 
A -1 (A -1 ViVj.F m n) and A~ 3 VjVjVfc-F mn . (ii) We trans- 
form the super-potentials into simple Poisson potentials 
in order to get rid of all derivatives acting directly on 
the sources, (iii) We insert the resulting quantities into 
the transverse traceless tensor A -1 ^™^, and perform 
some additional manipulations that lead to the final ex- 
pression. These steps are performed in detail in the fol- 
lowing. 

It is straightforward to expand the operator 7^ de- 



fined in Eq. (|A.6|I and to apply it to the source F k i given 
in Eq. (|27() . Taking into account the symmetry in the two 
indices k and I, inverting the tensor Poisson equation l)26|l 
then yields 

. TT _ a -1 (~TTkl jp \ 

n i:j -A (jy b kl ) 

= A" 1 ^ - p.-A- 1 (j kl F kl ) 

- 2A- 2 (^VfcV^i) + iA-^V, {i kl F kl ) 



At the lowest approximation, we may replace a by (—500+ 2^ 1 ^k^i (7 7 F, 

n ni\1 /O. 1 tt I 9. , mtl I 4\ t ,i . l . l . .11 . ... 



fafr) 1 ' 2 = 1 - U/c 2 + C(l/c 4 ). In the end, the elliptic 
equation for the lapse in the presence of hj^ is modified 
to 



-A-^iVjVfcVj (j mk r l F 



(B.l) 



\ Jhl7=o c 



TT 



(A.21) 



This is the desired CFC+ metric equation for the lapse 
function a. It can easily be transformed into Eq. (|29|l in 
Sect. 12.41 remembering that </> coincides with the CFC 
conformal factor at this level. Since U/c 2 is Newtonian, 
7 lfe 7^/!,T. T VfcVi?7/c 2 corresponds to the second post- 
Newtonian order for the dynamics. 

The equation for the shift can be obtained in princi- 
ple by contracting the 3-metric evolution with the help of 
the Euclidean metric 7 y . However, the new terms, by ref- 
erence to the conformally flat case, are proportional to a 
product of the type K (or (3) times hj^ '. They are therefore 



As the Poisson integral A~ 1 F k i converges, all other 
(super-)potentials entering Eq. I|B.1(I are also well defined. 
However, they cannot be easily handled. For instance, 
quantities such as A~~ 2 F mn or A F mn are a priori mean- 
ingless, which shows that the derivatives cannot commute 
with the integrals. In order to operate on the sources with- 
out meeting any serious restrictions, it is convenient to 
resort to the tool of Hadamard finite part regularization. 

It is out of our scope to review all properties of the 
Hadamard finite part. However, as it will be used inten- 
sively in the following, we recall its definition for com- 
pleteness, as well as those of its features we need in our 
derivation. When a function f(x ), x G R 3 , is smooth out- 
side a finite number of singularities, locally integrable, but 
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not integrable on K 3 , we can consider instead the new in- 
tegrand (\x — xo\/ro) B f(x), where B is a complex num- 
ber, r a positive number, and where \x — x \ denotes the 
Euclidean norm of x — x , Xq being an arbitrary vector 
of M 3 . The integral j\ x _ Xol>ro d 3 x (\x - x \/r Q ) B f(x) de- 
fined by means of the natural measure d 3 x = dx 1 dx 2 dx 3 
in Cartesian coordinates converges for B belonging to an 
appropriate domain D of the complex plane. It can be re- 
garded as a holomorphic function on D. It is then possible 
to extend If{B) = J d 3 x (\x — xo\/ro) B f{x) by analytic 
continuation as close to the point B = as desired, and 
to obtain its Laurent expansion T\^Ll f)k^ there, as 
explained bv lBlanchet fc Damourl l|l986|) . The zeroth or- 
der coefficient (i/)o is often referred to as the finite part 
integral of f(x), which in our notation reads 



FP B=0 



J d 3 x^f 



\X - Xq\ 

ro 



/(*) = (^)o. 



(B.2) 



The finite part integral of / may depend on the arbitrary 
radius ro; in fact, this will typically happen when the re- 
sult contains logarithms. Nonetheless, as long as / is inte- 
grable, its finite part integral coincides with J d 3 x f and 
does not show any dependence on r . 

An important property of the Hadamard regulariza- 
tion is that the finite part Poisson integral of a smooth 
function always exists, whereas the simple Poisson inte- 
gral may not. The covariant expression of this finite part 
reads 



FPA- o 1 / = FP B=0 



dVV7 f\x'~x \\ B f(x') 



-4n 



ro 



(B.3) 



where the Euclidean volume element has been written as 
d a3-\/7 in an arbitrary coordinate system. When A -1 / 
exists, it satisfies FP A^ 1 / = A -1 /. In any case, the reg- 
ularized Poisson integral is the particular solution of a 
Poisson equation of the type Ag = /. Thus, the operator 
FP A -1 constitutes a genuine generalization of the ordi- 
nary Poisson operator A -1 and will be denoted as A" 1 
henceforth. It has the important property that it com- 
mutes with the spatial derivatives Vj, which allows us to 
work on the form of the elementary (super-)potentials by 
applying simple and systematic rules. 

An important formula related to the Hadamard regu- 
larization is the one giving the generalized Poisson integral 
of the distance to the field point x, 



\x — x 



i\A 



\x — X 



l\A+2 



(A + 2){A + 3Y 



(B.4) 



for an arbitrary complex exponent A =/= —2, —3. Notably, 
it can be used to evaluate the action of the operator A^ 1 = 
FP c =o/ dV y/j/(-±nr% \x-x"\) on the "r Q " -potential 
FP B =o/dVV7|a; - x'\ a+B f /(-4irr B ) with a e Z. By 



permuting the two triple integrals, we obtain the relation 



A" 1 FP B=0 



= FP s =o ■ 



„na+B+2 



f 



-4nr B (a + B + 2) (a + B + 3) 



(B.5) 



The result has the same form as the source. If we make 
A" 1 act on it, we arrive at a quantity of the same type. 
This provides a straightforward procedure to determine 
the action of A~ p = (A~ 1 ) p , p 6 N, on the original inte- 
gral iteratively. In this way, we find: 



A- p FP B=0 



FP B=0 



d 3 x'Vj 
-4nr B 



J\a+B 



f 



(a + B + l + ip)- 1 



x (a + B + 2pY x ...(a + B + 2y 



X IX — X 



l\a + B+2p 



(B.6) 



If the usual "r°" -potential with source / exists, the pth it- 
erated Poisson integral A" 1 of its pth derivative also does. 
It is equal to V ^ . . . V, applied to the right-hand side of 
Eq. ljB.6|l . When we put all derivatives under the integra- 
tion symbol, we end up with a convergent integral. At this 
stage, one does not need any regularization anymore, so 
we may take B = if none of the terms (a + 2p + 1), 
(a + 2p), . . . , (a + 2) vanishes. Finally, we pull the coef- 
ficients (a + 1 + 2p)- 1 {a + 2p)~ x . . . (a + 1 + p)" 1 out of 
the integration symbol and reintroduce the finite part. We 
have then proved the formula giving the explicit action of 
A~ p on arbitrary sources, 



d 3 x'^j ] 

X 

V n . . . V lp FP S=0 / ~ xT +B+2p f 



(a + l + 2p)(a + 2p)...(a + 2) 



(B.7) 



to be valid for a + 2p + 1 ^ N. 

This allows us to write the potentials and super- 
potentials in a form suitable for numerical integration. 
For this purpose, we can recourse to fairly standard tech- 
niques, some o f whic h have been used in particular by 
Blan chet et al to deal with the derivative of the 

Newtonian super-potential. With the help of Eq. l|B.7(l 



and the commutation relation A^ Vj = Vj A x , we trans- 
form the expression of the (super-)potentials A - VjVji^fci 
or A ;l V,V , V/ V;/-„,„ entering the non-conformal part of 
the 3-metric (with possible index contractions). By doing 
this, we try to minimize the number of free indices that re- 
main inside the integral and, for numerical reasons, to get 
rid of all spatial derivatives of the densities D* or S* in the 
final expressions. Let us detail for instance the transfor- 
mation of K~ 2 V tV jFki- It is achieved in four steps: (i) We 
let both derivatives commute with the integration symbol, 
A- 2 V l V JJ F 1 fci = ViVjA~ 2 F k i, (ii) we rewrite A~ 2 F ki ac- 
cording to Eq. (|B.6|) with p — 2, (iii) we make one of the 
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derivatives act on the kernel \x — x'\ 1+B so that there re- 
mains only one unevaluated spatial derivative operating 
on a linear combination of simple Poisson integrals, and 
(iv) we let the derivative act. This yields 



A 2 ViVjF k i = - 



la 



kl + IjmX 



ViA-^j^Fu) 



(B. 



The respective transformation of A ViVjVfc-F mn is very 
similar: 

A 3 ViV : ,Vfei ;i mn 

= ^V i V J -V fe FP B=0 f ^/g.\ x -xf+ B F mn 
24 ,/ -47rrn 



1~ (dVV7 

8 Vfe ./— 



{x p -x' p )(x q -x' q ) 



lipljq- 

+ %\x - x' 



x — x' 

TP 

- 1 - mil 



&Y(ij%)pri l 'A a , l F mn - 37 (ij 7 fe )pA :i , 1 (x p F mn ) 
-2%^ j)qX Pij k A x 1 {x q F mn ) 

+ j ip j 3q VkA~ 1 {x p x q F mn ) 



(B.9) 



We are now able to deduce the expression for up to 
second post-Newtonian order, parametrized by means of 
the intermediate elementary potentials S, Si, %, IZi, and 
Sij as 



hJ7 = Is 



- 3x fc V (j 5 j)fc + jW m Vi (r l Ski) 
-x k x l ViVjS k i + 3\7 {l S j} - -x^iVjSk 



Hi 



\l kl S kl + x k i m \J m S kl - 7 fei V fe 5 ; 



(B.10) 



This is the expression we use in Eq. (|35[1 in Sect. 12.41 with 
the potentials determined by the scalar/vector/tensor 
Poisson equations 



Appendix C: Multipole expansion of the 
intermediate potentials 

Here we derive the multipole expansion of the interme- 
diate potentials S, Si, %, IZi, and Sij needed for ob- 
taining boundary conditions for the se potentials. This is 
done b y means of the formula (C.9) of B lanchet fc Pouiadel 
(2002) speci alizing the matching relation first established 
in iBlanchetl l)l998|) for retarded quantities. 

For any generic source / which admits outside the 
system a multipole expansion of the form M(f) = 



fp( n ) rP w ith n = x/r and p n < —2, the mul- 
tipole expansion of the Poisson integral A -1 / is given by 



M 



(A-V)-t!^ v„(I) 

+ A 1 X(/). 



x n 'f 



(C.l) 



In the special case where the source / has compact sup- 
port, M. (/) is identically zero and thus the last term above 
vanishes. We recover the standard multipole formula used 
in electrostatics for spatially limited systems. 

At this stage, the multipole expansions of all our el- 
ementary potentials may be derived by application of 
Eq. IjC.ll) . We start with S^ which goes to zero at the 
highest order. It involves the monopole integral with 
non-compact-supported source J d 3 x iU\7 jU / (— Air). 
Remarkably, this integral admits an alternative expression 
whose source has compact support, which is a very use- 
ful feature for the numerical calculations. To perform the 
transformation, we first replace the second potential U in 
the integrand V i {7V J {7/(— 47r) by Ax/2, where \ denotes 
the Newtonian super-potential \ — Jd 3 a;'|a; — x'\D* . 
Next, we integrate this Laplacian by parts, being care- 
ful of keeping all contributions from the derivatives of r B . 
The result is 

FP S=0 / ^5r*V, ; 



-Airr B 



(!) 



V,-A£7 



FP B=C 



d 3 X y/j 



-4:irr B 



B(B + l)r B - z VjU 



2Br L 



(!) 



(C.2) 



We remark here that the second finite part on the right- 
hand side vanishes. Indeed, the integration does not gen- 
erate any pole in B able to compensate the cancellation 
of the pre-factor B itself. Consequently, after a last inte- 
gration by parts, we arrive at the equality 



d 3 x V7 

-47T 



ViUVjU 



1 



(C.3) 



It is not difficult to check that VyX = 
7 ifc x fc A- 1 (-47rD*) - A~ 1 (-4:%j jk x k D*) by virtue 
of relation i|B.7|) . Using an integration by parts of the 
type 



J d 3 x ^fA~ 1 g = J d 3 x^gA- 1 
we are finally able to show that 

-4tt 3 



f, 



(C.4) 



= - I d 3 X y/^D* [ ^kX k V t U+-%U 



(C.5) 
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Inserting the latter relation into the general formula i|C.l(l 
specialized for / = -4ttS*S*/D* - V.UVjU, we get 



M{S l3 ) = - 
r 



M -4tt- 



d 3 ai sf^D* 
S*S* 



D* 



D* 



When examining the second term in the above equa- 
tion, one immediately observes that M. (S* S* / D*) = 
as the source (S*S* /D* 2 )D* has compact support. 
Furthermore we notice that M(V iUV jU) = 0(l/r 4 ), 
hence Aq 1 M(V i UW j U) = 0(l/r 2 ) according to dimen- 
sional analysis. It may also be checked directly with 
the help of the Matthew formula Aq (n <l1 



-,ii> r a\ _ 



,<il 



^>r a+2 /[(a + 2-l)(a + 3 + l)]. 



We continue with the calculation of M.(Si). The part of 
its source with non-compact support amounts essentially 
to x k \7 iUV k U . The integral of the latter quantity may be 
transformed with the same technique as the one used to 
establish Eq. QClfy . We find that 



FP_B = o 



• FP B=0 



d 3 x V7 



+ 2 7 ^V fcX V / V l C/V (C.7) 



We may then perform an integration by parts affecting 
the derivative ViD* of the first term and the derivative 
VfeX °f the second term. The contribution of V^r B is pro- 
portional to B. It can give rise to a definite non-zero re- 
sult only when it is multiplied by the factor 1/B coming 
from the radial integration of 1/r 3 . As the correspond- 
ing angular integral vanishes, so does it. Noticing that 
AxVjt/ = 2U\7iU — \?iU 2 , we see that the resulting in- 
tegral reads 



1 



jD* (x k 



V 4 V fc x + v,x 



- FP S=0 / d 6 x ( — 



V ? ;C/ 2 . 



(C.8) 



After integrating the second term by parts, we get an in- 
tegral whose source is proportional to B. Following the 
same argument as before, it can be shown to be simply 
zero. Expressing again the super-potential as a combina- 
tion of Poisson integrals, and performing other integra- 
tions by parts including those of the form (|C.4Jl , we arrive 
at 



FP B=0 



= - / d A x 



jD*% (x j U + x j x k V k U 



(C.9) 



On the other hand, the source i V^V^ behaves as 
M 2 ^ ik n k /r 3 = -M 2 Vi{r- 2 /2), where M is the bary- 
onic rest mass density defined in Sect. 12.41 A direct inte- 
gration leads to the relation AQ 1 (l/r 2 ) = lnr + const., 



which can also be seen from a dimensional argument, 
so that Aq 1 M{x k V iUV k U) decreases asymptotically as 
— M 2 Vi(lnr)/2. Therefore, we have proved the approxi- 
mate equality 



(C.6) M{S l 



1 



x k S i S k +%x j (U + x k V k U) 



M 2 , /n ( 1 
2r \ r z 



(CIO) 



In the multipole expansion of 7^, the only term en- 
tering its composition that is generated by a quan- 
tity with non-compact support is the Poisson potential 
A -1 (lijl kl x^V ' k UV 'iU) . Its monopole part is composed 
of an integral over the source, which simplifies by virtue 
of the relation 2j kl V k UViU = AU 2 - 2UAU to 



FP, ; -u / {ji^x^VkUVtU 



J d 3 Xy/^D*J lk X k U, 



(C.ll) 



and of a contribution originating from the Poisson 
inverse operator applied to M^ij^x 3 V k UViU) = 
-M 2 V;(r- 2 /2) + 0(l/r 4 ). We have thus 



M 2 



Ao 1 A<(7ii7 fcl ^VfcJ7VjJ7) = - — 



2r 



(C.12) 



The treatment of the leading term in the multipole ex- 
pansion of IZi is similar. The integral of V i(x k x l V k UV ill) 
is simply zero, by virtue of the Gauss theorem, and the 
monopole part of — A -1 ^ i{x k x l V k UV iU) is the same as 
that of 2^ kl x : >V k UV l U. 

We conclude with the multipole expansion for the po- 
tential S. Since its source is already of compact support, 
it can simply be expanded with the help of the standard 
multipole formula mentioned after Eq. IjCljl . 

In the end this yields the desired multipole expansions 
of the elementary potentials in the second second post- 
Newtonian expansion, which correspond to Eqs. l|36H40fl 
in Sect. 



Appendix D: Post-Newtonian and 

post-Minkowskian link for gravitational waves 

The equation for null outgoing geodesies in the ADM 
gauge reads t-r/c - 2A/ A DMGln[r/(c?7)]/c 3 + 0(G 2 ) = 
const., with M A dm = M + 0(l/c 2 ) being the ADM mass 
and J] an arbitrary positive constant with the dimension 
of time. For r large enough compared to the typical wave- 
length, this also gives the link between the time t and the 
distance r of a field point x located near I + . In the post- 
Newtonian framework, r/c is regarded as a small quan- 
tity (with respect to t), which precisely reflects the fact 
that the post-Newtonian metric is not accurate in the 
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wave zone. Therefore, the expression (|35[1 for hj? can- 
not be used there. Instead, we should consider the post- 



Minkowskian version of Eq. I|A.15|) . where all quantities 
are expanded in powers of the gravitational constant G. 
Higher order terms in 1/c 2 must not be a priori neglected 
as long as they appear at a level of approximation below 
the one we want to reach. 

In general, this approach is not guaranteed to lead to 
a well-defined perturbative scheme in the ADM gauge in 
the sense that the contribution of order G n+1 may become 
greater than the contribution of order G n in the neighbor- 
hood of I + . It is worth to briefly point at the origin of the 
possible problem. Once the field hnj^ is known at the 
nth post-Minkowskian (nPM) level, the next approxima- 
tion h n +ij^ is determined by the wave-like equation of 



the type Uh n+1 ^ = ^A n+lkl (h^)+0(G n + 2 ) ob- 
taincd by combing Eqs. (|A.4ll and (|A.5(1 as explained in 
Appendix lAl The source A n+1 . . depends non-linearly on 



TT 



hlj , which can thus be replaced by its nPM value h nkl . 
In the end, the solution of the previous wave equation 
is the retarded integral of j]j Tkl A n+ i kl (h TT ), which we 

denote in short as □ _1 7?™A n _|_i fei . The transverse trace- 
less projector can itself be written as a converging integral 
with the help of a modified version of Eq. I|B.7(I , in which 
the derivatives are applied to the kernel |x — x '\ a+B+2p 
with a = —1. Even if the retarded integral converges, 
it may be conveniently replaced by its Hadamard finite 
part. The numerical result is not affected, but the inte- 
gral appearing under the operator FPs=o can now be 
shown to commute with the integral defining h n+ ij^~ after 
an appropriate analytic continuation. As a consequence, 
the finite part of 7y T regarded as an integral operator, 
which may be referred to as 7^™ without any ambiguity, 
can be pulled out from the d'Alembertian inverse sym- 
bol: D-^™A n+lkl = j™D^A n+lkl . The behavior 
of A„ + i i:) (/i„ TT ) when r — » +00 follows from a formula 
extending Eq. (|C.1(1 to the case where the quantity to be 
calculated is a retarded integral. In particular, the multi- 
pole expansion of hj^ involves the retarded integral of the 
multipole expansion of the source, regularized by means 
of the Hadamard finite part Dq 1 A4(A n+ i i j). This contri- 
bution turns out to generate terms tending towards zero 
as (In r) 2 /r at the 2PM order, which are thus out of con- 
trol. The structure of the field has not been investigated 
yet at higher PM order in ADM gauge, but a similar phe- 
nomenon is likely to happen in that case too. 

Fortunately, the previous problem can be cured at the 
post-linear level either by moving to some radiative gauge 
where u = t — r/c is a conformal coordinate, or by ab- 
sorbing the logarithms into the arguments of some multi- 
pole moments. In harmonic gauge, the logarithms occur- 
ring in the post -Minkowskian iteration have been shown 
l)BlanchetHl987h to be removable to all orders. In ADM 
gauge, the possibility of such an eliminat ion has been 
checked explicitly at the post-linear level bv lSchafe'il l)l990|) 
from a formula adapting the derivation of Appendix lAl in 
order to include in the Hamiltonian the relative Newtonian 



TT 



part of the 2PM corrections to all terms quadratic in h. 
(or 71"^). The new wave equation for the field variable 
reads 



ah TT = -i™ 



■T n V m UV n h kl ~ru^ kl 




11 ' UAh TT 
0((^ TT ) 2 ).(D.l) 



The Green function method provides the required solu- 
tion, but the computation is rather delicate. The result 
is a retarded integral, on which it remains to apply the 
transverse trace- free projector 7y . The action of the 

second spatial derivatives V^V; on a quantity of the form 
f(t—r/c)/r leads to the "monopole" term A fkpliqri p n q f (t— 
r/c)/(rc 2 ), each dot denoting a time derivative, plus 1/r 2 
corrections. It follows that A _1 VfeV;[/(t — r/c)/r] — 
lkpli q n p n q f(t — r/c)/r + 0(l/r 2 ). Similarly, it is not 
difficult to check that A~ 2 V ;V 3 -V fcV ;[/ '(t - r/c)/r] = 
lipljglkrli s n p n q n r n s f (t - r/c)/r + 0(l/r 2 ). This entails 
that the operator (|A.6|I applied to a term f(t — r/c)/r re- 



duces to the projector P i 



TTfci 



Using the latter property, 



we recover the well-known relation 



1 



-P. 



x <M M 



TT kl 



r 

t 

c 



2G 
IT 



—t- \ M + O 



-0(G 2 



•0 4 



c 2 ) ) ln (ct7, 



(D.2) 



with Mij(t) = J d 3 x y/jFij/(— 4n) being the Newtonian 
monopole moment associated to the source Fij . 

Because of the algebraic transverse traceless projec- 
tion, we may substitute M < ij > by Mij without affecting 
the outcome in the above relation. Since we still work at 
the 2PN level, we do not need to include any higher order 
corrections to the expression of Mij. The retarded argu- 
ment t-r/c- 2GM ADM /c 3 \n[r/(c7j)} + 0{G 2 ) coincides 
with the null radiative coordinate. It can be viewed as the 
actual retarded time i rot of the field variable at I + up 
to the 1PM order. Taking the preceding into account, we 
arrive at 



,TT 



R 



TT/cZ 



-^<fei>(iret) 



O 



1 



(D.3) 



We now improve the formula for M^ by making use of 
Eq. EH) . This leads to 

(D.4) 

Modulo 1PN errors, this integral equals 
4G / d 3 z y/iD(%ai j v <k v l> + x k % <t V 3> U)/c A . This 
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is proportional to (twice) the second time derivative of 
the Newtonian quadrupole 2Q t j at this level, as can be 
checked with the help of the relation 



^ J d 3 x^D*f(x,t) 

d 3 x^D* UWi + dA /(*,*), 



(D.5) 



which is valid for any C 1 function /. We have thus derived 
the Newtonian quadrupole formula 



lTT pTTfci 



2GQki 



(D.6) 



in the neighborhood of future null infinity. 

It now is possible to investigate the asymptotic behav- 
ior of hfj^ in a similar way as we did for the intermediate 
potentials. The computation does not present any partic- 
ular subtleties in contrast to that of the full field variable 
dominant contribution. The operator commutes 

with the Poisson integral so that Al(A~ 1 7^ Tfc ' Fu) — 
Z/J^Mik^Fki) = 7™(M fc i/r) +0(l/r 2 ). The action 
of the transverse trace-free projector is computed with the 
help of Eq. (|B.4|) . This gives 



3 7fo-7 fe ° + jS^ j)p n l n p 



1G 



Y^ A /i P l 3q n p n q n k n l 



5 

~ 16 

M k i(t) 



(%jn k n l + 'j kl % A / jq n p n' 1 ) 



(D.7) 



2GQij/c A . By applying the operator P^, 



jTTfcii^PN 



with M. 

on both sides of Eq. ijD.7(l . we find that n kl 
P^j™GQki/ (4rc 4 ). As a conclusion, with the notations of 
Sect.E3 we have (cf. Eqs. J52E2)): 



jTTfcii 2PN 

ij n kl 



(D. 
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